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Abstract 

A  computational  aeroelastic  analysis  of  a  micro  air  vehicle  (MAV)  is  conducted.  This 
MAV  has  a  24  inch  wing  span,  and  is  designed  for  local  area  reconnaissance.  Wind  tunnel 
data  for  the  MAV  with  a  rigid  carbon  fiber  wing  and  a  flexible  carbon  fiber  ribbed  nylon 
wing  are  compared  to  CFD  results  incorporating  static  and  dynamic  deformations.  We  use 
laser  vibrometry  to  determine  the  mode  shapes  of  the  flexible  wing.  From  these  shapes,  CFD 
grid  deformations  are  calculated  as  part  of  a  closely-coupled  aeroelastic  solution  method. 
The  accuracy  of  MAV  performance  predictions  using  CFD  with  and  without  aeroelastic 
modeling  is  evaluated  against  previous  wind-tunnel  experiments.  The  performance  benefits 
of  the  flexible  wing,  and  the  applicability  and  limitations  of  the  model  are  evaluated  in  the 
present  research  effort.  Some  suggestions  are  made  as  to  improvements  that  can  be  made 
to  increase  the  range  of  applicability  and  the  accuracy  of  the  model. 
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Computational  Aeroelastic  Analysis  of  Micro  Air  Vehicle  with 
Experimentally  Determined  Modes 

I.  Introduction 

Micro  air  vehicles  (MAVs)  are  small  unmanned  aircraft  developed  for  a  variety  of  spe¬ 
cialized  missions  in  Intelligence,  Surveillance,  and  Reconnaissance  (ISR).  According 
to  the  Defense  Advanced  Research  Projects  Agency  (DARPA)  the  MAV  is  a  small  aerial 
robot  capable  of  exploiting  its  small  size  and  mobility  for  a  broad  range  of  remote  operations 
(1141).  The  MAV  must  provide  a  stable,  inexpensive,  and  expendable  surveillance  platform 
(T2T|) .  Similarly  to  the  gun  platforms  provided  by  larger  aircraft,  the  MAV  must  provide 
precise  directional  control  for  mission  success,  and  it  must  do  this  in  a  very  demanding 
flight  regime.  Because  they  have  small  dimensions  and  must  fly  at  low  speeds,  the  Reynolds 
numbers  experienced  by  typical  MAV’s  range  from  104  to  105.  At  such  low  Reynolds  num¬ 
bers  the  vehicle  is  susceptible  to  sudden  separation  and  loss  of  lift  due  to  very  large  gusts. 
On  a  typical  day  the  wind  may  vary  by  10  mph,  which  is  a  significant  velocity  component 
for  an  aircraft  traveling  at  10  to  30  mph  (]10j). 

MAVs  are  highly  energy  and  weight  constrained.  Because  of  this,  all  components 
of  the  MAV  must  be  more  tightly  integrated  than  those  of  a  traditional  aircraft,  and  ex¬ 
ploiting  multi-functionalities  across  components  becomes  crucial  to  achieving  useful  MAV 
performance  (1I4|).  Of  particular  interest  along  these  lines  is  the  impact  on  performance  of 
a  highly  flexible  wing  at  low  Reynolds  number. 

The  flexible  wing  is  desirable  for  a  variety  of  reasons.  Logistically,  the  flexible-wing  has 
a  smaller  packaging  footprint  which  is  crucial  for  its  success  as  a  man-portable  platform. 
Aerodynamically,  it  passively  adapts  to  local  flow  disturbances,  providing  more  constant 
lift  in  gusty  conditions.  This  process  is  termed  “adaptive  washout”,  and  acts  by  the  gust 
producing  a  deformation  in  the  wing  which  reduces  its  lifting  efficiency,  but  the  higher 
dynamic  pressure  in  the  gust  leads  to  near  constant  lift  production  (TIP).  The  dynamic 
response  of  the  wing  structure  must  be  tuned  to  provide  the  correct  amount  of  deformation 
with  gust  onset,  while  at  the  same  time  meeting  packaging  and  strength  requirements. 
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The  already  recognized  importance  of  fluid-structure  interaction  in  large  aircraft  de¬ 
sign  becomes  a  dominant  consideration  for  the  MAV,  which  is  especially  sensitive  because 
of  its  light  weight  and  large  surface  to  volume  ratio  (1141).  The  flexible  wing  thus  provides 
a  challenging  opportunity  to  exploit  aerodynamic,  logistic  and  control  benefits  in  the  MAV 
design. 

1.1  Research  Objectives 

The  goal  of  this  research  is  to  implement  a  closely-coupled  aeroelastic  model  using  a 
combination  of  computational  and  experimental  tools.  The  model  will  be  evaluated  against 
a  wind  tunnel  characterization  of  a  rigid  and  flexible  wing  MAV.  The  structural  model  will 
rely  on  experimentally  determined  mode  shapes  which  are  coupled  with  a  time-accurate 
CFD  solution.  The  model  addressed  in  this  paper  allows  the  unsteady  behavior  of  the  MAV 
to  be  evaluated,  and  the  quasi-static  deformation’s  effect  on  performance  can  be  compared 
to  that  measured  in  the  wind-tunnel  testing. 

1.2  Prior  Experimental  Work 

The  MAV  configuration  investigated  in  this  thesis  was  evaluated  in  previous  wind- 
tunnel  tests  by  Deluca  (7J.  The  flexible  wing  MAV,  shown  in  Figure  |JJ  has  a  24  in  wing¬ 
span,  and  is  constructed  of  carbon  fiber  and  parachute  nylon.  The  wing  has  a  carbon  fiber 
leading  edge  and  chord-wise  carbon  fiber  ribs  spanned  by  the  parachute  nylon,  creating 
some  resemblance  to  a  bat  wing.  The  rigid  wing  MAV  shown  in  Figure  [T]  is  similar  to  the 
flexible  wing  version  except  that  the  wing  is  constructed  entirely  of  carbon  fiber  (further 
physical  details  in  Table  HJ. 

The  results  presented  in  (7J)  are  time-averaged  force  measurements,  and  establish  the 
static  stability  of  the  aircraft.  The  wind  tunnel  tests  show  that  the  static  deformation 
of  the  flexible  wing  at  high  angle  of  attack  increases  the  maximum  lift  to  drag  ratio,  as 
well  as  maximum  lift  coefficient  over  the  rigid  wing.  Also,  as  noted  by  DeLuca  (115j).  the 
unsteady  behavior  of  a  possible  laminar  separation  bubble  just  aft  of  the  peak  suction  has 
a  significant  affect  on  the  MAV’s  performance  and  is  very  sensitive  to  flow  disturbances  or 
small  vibrations  of  the  wing. 
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Figure  1:  Flexible  and  rigid  wing  MAVs 


Table  1:  MAY  physical  properties 


Flexible  wing 

Rigid  Wing 

Mass 

0.320 

0.360 

kg 

Plan- for  nr  Area 

6.032  x  10-2 

6.032  x  10-2 

2 

m 

Root  Chord 

0.1524 

0.1524 

m 

Mean  Aerodynamic  Chord 

0.1067 

0.1067 

m 

Span 

0.6096 

0.6096 

m 

Leading  edge  thickness 

6.350  x  10"4 

6.350  x  10"4 

m 

Membrane  thickness 

1.270  x  10”4 

1.270  x  10"4 

m 

Aspect  Ratio 

6.1604 

6.1604 
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Other  researchers  have  conducted  wind  tunnel  testing  of  MAVs  with  membrane  wings 
which  have  shown  similar  improvements  in  performance  with  the  wind  tunnel  results  pre¬ 
sented  in  (H).  A  slightly  smaller  series  of  vehicles  examined  in  (I21j)  show  the  benefit  of 
membrane  over  rigid  wings  in  stall  characteristics.  Another  interesting  aspect  of  membrane 
wing  performance  noted  in  (1211  is  the  adaptive  washout  of  the  flexible  wing  during  gust  and 
at  high  angles  of  attack.  This  provides  improved  performance  by  smoothing  the  production 
of  lift,  while  at  the  same  time  decreasing  the  weight  of  the  aircraft  over  that  of  the  rigid 
wing. 

1.3  Computational  Modeling  of  Membrane  Wings 

A  series  of  slightly  smaller  rigid  and  flexible  wing  MAV’s  are  investigated  in  (111!)  using 
a  closely-coupled  method  incorporating  a  Pressure  Implicit  with  Splitting  Operators  (PISO) 
algorithm  to  solve  the  incompressible  Navier-Stokes  equations  and  a  hyperelastic  Moody- 
Rivlin  model,  based  on  the  approach  presented  in  (T12T) .  They  have  shown  that  the  billowing 
of  the  membrane  can  have  an  offsetting  affect  on  lift  production  to  the  adaptive  washout,  so 
the  lift  might  be  little  changed  from  the  rigid  wing.  Their  results  also  show  high  frequency 
membrane  vibration  in  steady  free-stream  (TllTt.  The  membrane  deformations  at  nominal  to 
high  angle  of  attack  tend  to  increase  the  maximum  lift-to-drag  ratio. 

1-4  Research  Approach 

The  implementation  and  validation  of  the  aeroelastic  model  begins  with  steady,  non¬ 
deforming  solutions  to  establish  a  base-line  comparison  with  the  rigid  wing  wind-tunnel 
results.  Static  structural  deformations  are  compared  to  the  flexible  wing  wind-tunnel  mea¬ 
surements  to  evaluate  the  validity  of  the  linear  modal  model.  The  flexible  wing’s  effect  on 
dynamic  performance  is  measured  by  a  simple  gust  response. 

Particular  attention  is  paid  to  the  contributions  to  lift  variation  of  the  flexible  wing 
under  steady  flight  conditions  as  well  as  gust  since  lift  production  of  the  MAV  is  a  critical 
factor  in  its  performance  and  mission  accomplishment. 
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1.5  Thesis  Organization 

Chapter  fU introduces  the  MAV  concept,  prior  work  in  the  field,  and  the  current  re¬ 
search  approach.  Chapter  jXJ]  presents  a  detailed  description  of  the  solution  methods  and 
the  organization  of  the  overall  solution  algorithm.  The  results  are  presented  in  Chapter  III. 
with  conclusions  following  in  Chapter  lIYl 
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II.  Research  Methods 


In  general  the  computational  aeroelastic  problem  is  a  three  field  problem,  with  conserva¬ 
tion  equations  governing  the  fluid  and  structure,  and  some  arbitrary  choice  of  governing 
equations  for  the  mesh  geometry  in  each  zone  of  the  simulation  (T12).  There  are  several 
approaches  to  this  type  of  problem.  In  a  closely-coupled  solution  the  surface  forces  are 
transfered  between  the  fluid  solution  and  the  structural  model  at  each  time  step  or  itera¬ 
tion,  while  a  loosely- coupled  solution  performs  updates  of  deformations  at  a  less  frequent 
interval  (fl3l) .  In  a  fully-coupled  solution  the  equations  of  motion  for  the  fluid,  structure, 
and  grid  are  reformulated  into  one  set  of  governing  equations  (0]).  Figure  E)  shows  the  or¬ 
ganization  of  the  high  level  solver  components  in  a  computational  aeroelastic  analysis.  The 
current  method  is  closely-coupled ;  it  communicates  loads  between  the  fluid  and  structure 
and  updates  the  grid  at  each  time  step. 

No  matter  the  solution  approach,  communication  between  the  computational  domains 
of  the  fluid  and  structural  simulations  is  a  limiting  factor  on  the  accuracy  of  the  aeroelastic 
analysis.  The  problem  of  communication  between  the  problem  domains  can  be  eliminated 
by  making  the  fluid  and  structure  meshes  point-matching  along  the  surface  interface.  This 
requires  the  preliminary  analysis  of  the  problems  to  be  very  closely  coordinated  and  forces 
somewhat  unreasonable  constraints  on  the  two  grids.  The  computational  meshes  required 
to  accurately  and  efficiently  solve  the  two  problems  differ  greatly  in  complexity.  CFD  grids 
are  usually  stationary,  and  the  equations  of  fluid  dynamics  are  posed  in  an  Eulerian  frame. 
The  CFD  meshes  are  also  highly  stretched  near  solid  surfaces  to  capture  gradients  due  to 
viscosity.  The  structural  equations  of  motion  are  usually  posed  in  a  Lagrangian  reference 


j  Spline 

Spline 


Structure 


Displacements 


Figure  2:  General  aeroelastic  algorithm 
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frame,  and  the  grid  deforms  with  the  structure.  This  mesh  tends  to  be  much  more  isotropic 
than  the  CFD  mesh,  and  the  resolution  requirements  for  acceptable  accuracy  are  much  less 
stringent.  This  is  especially  true  in  the  case  of  a  modal  analysis,  where  a  very  simplified 
model  of  a  complex  aircraft  geometry  is  often  sufficient  to  accurately  represent  the  lower, 
energy  containing,  modes  of  dynamic  motion  (QQ. 

The  structure’s  time  response  can  be  calculated  by  directly  integrating  the  full  equa¬ 
tions  of  motion,  or  integrating  the  equations  on  a  reduced  basis  of  the  lower  modes.  The 
modal  approach  has  the  advantage  of  requiring  only  one  solution  of  the  structural  prob¬ 
lem,  and  the  flexibility  of  using  either  computationally  or  experimentally  determined  mode 
shapes  and  frequencies.  Alternatively,  the  direct  integration  approach  requires  the  full  struc¬ 
tural  simulation  to  be  solved  concurrently  with  the  fluid  dynamics  calculations.  This  added 
cost  is  hardly  justified  for  problems  dominated  by  “slow”  physics,  that  is,  the  dominant 
frequencies  of  the  structural  response  are  much  lower  than  the  characteristic  frequencies 
of  the  flow.  Thus  a  modal  analysis  takes  advantage  of  two  efficiencies.  The  mesh  size  is 
reduced,  and  only  one  solution  or  experiment  for  the  structure  portion  of  the  aeroelastic 
problem  is  needed. 


2.1  Fluid  Solution:  Navier-Stokes 

The  Navier-Stokes  equations  are  a  continuum  formulation  for  the  conservation  of  mass, 
momentum,  and  energy  in  a  fluid.  The  integral  form  of  the  governing  equations  is 


d_ 

dt 


W  dV  + 


(J7  —  1FV)  ■  dA  =  0 


(1) 


Where  the  vector  of  conserved  variables,  and  the  fluxes  are 


p 

pU 

0 

pu 

pull  +  pnx 

^x^~xx  TlyTxy  ^z^~xz 

pv 

T  = 

pvU  +  pny 

Ev  = 

^x^~yx  H-  'fl'y^-yy  +  'fl'z'l'yz 

pw 

pwU  +  pnz 

'H'x'Tzx  'H'y'l'zy  ^z^~zz 

_  PE . 

(pE  +p)U 

Tlx®x  Tly®y  ^lz®z 
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The  contravariant  velocity,  U,  is 


U  =  v  ■  n  =  unx  +  vny  +  wnz 


(3) 


0  in  Equation  [2]  is  given  by 
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utxx  +  vrxy  +  wtxz  +  kdT/dx 

UTyx  +  VTyy  +  WTyz  +  kdT/dy 

utxx  +  vrxy  +  wtxz  +  kdT/dx 


The  shear  stress  tensor  (following  Stoke’s  hypothesis)  is  given  by 


'  xy  —  1  yx 
2 

~  31 


du 

dv 

dw 

dx 

dy 

dz 

f  du 

dv 

=  p( 

KdV 

dx 

dv 

du 

dw 

dy 

dx 

dz 

'  dw 

du 

=  p{ 

dx 

+  lh 

dw 

> 

dv 

du 

dz  dy  dx  J 
f  dv  dw\ 

TyZ  =TZy=y^—  +  —j 


(4) 


(5) 

(6) 

(7) 

(8) 

(9) 

(10) 


An  assumption  of  incompressibility  (constant  density)  leads  to  a  simplified  form  of  the  gov¬ 
erning  equations,  which  decouples  the  conservation  of  energy  from  the  mass  and  momentum 
equations. One  method  of  solving  this  set  is  elaborated  in  section  12.1.11 


2.1.1  Segregated  Incompressible.  The  low- velocity,  segregated,  semi-implicit  method 
for  pressure-linked  equations  (SIMPLE)  method  is  an  iterative  procedure  to  solve  the  in¬ 
compressible  Navier-Stokes  equations.  It  applies  a  pressure  correction  to  the  decoupled, 
incompressible  continuity  and  momentum  relations  (1 19).  The  pressure  is  written  as 


P  =  Po  +  P 


(11) 
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where  po  is  the  guessed  value  of  pressure  and  p  is  a  pressure  correction.  Similarly,  the 
velocity  components  are  written  as 


u 

Uo 

u 

V 

= 

Vo 

+ 

V 

w 

w0 

w 

(12) 


The  velocity  corrections  are  computed  using  an  approximate  form  of  the  momentum  equa¬ 


tion 


du  _  dp 

P  dt  ~  dx 


n 9v 

P  dt 

o— 

P  dt 


dp 

dy 

dp 

dz 


(13) 


The  velocity  correction  is  assumed  to  be  zero  at  the  previous  iteration  step  so  the  left  hand 
side  of  Equation  J3]  can  be  simplified  as 


un+1  -  un 


At 


=  P- 


,n+ 1 


At 


(14) 


Substituting  Equation  UJ]  into  HD  gives  the  current  velocity  corrections  as 


v  =  -A%  (15) 

»  =  ~A% 


where  A  is  the  pseudo-time  step  divided  by  the  density.  Substituting  [12]  and  CCS  into  the 
continuity  equation  gives 


duQ  dvo  dwo\  ( d2p  d2p  d2p\ 

foi  +  ^  +  ^)+A{e^  +  w  +  d^)  = 


(16) 


which  is  solved  for  the  pressure  corrections,  which  are  then  used  to  update  the  pressure  and 
velocity  field.  A  central  difference  is  generally  used  in  equation  [13]  for  the  p  derivatives, 
which  can  result  in  odd-even  decoupling,  for  this  reason  the  SIMPLE  method  is  generally 
solved  on  a  staggered  grid  to  prevent  odd-even  decoupling  in  the  pressure  field  (TI9]).  We 
use  the  commercial  CFD  package,  FLUENT,  which  implements  the  above  algorithm.  The 
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velocity  gradients  are  found  using  either  a  first  or  second  order  upwind  discretization  (in 
FLUENT).  Both  of  the  schemes  are  used  in  the  present  research.  The  effect  of  the  spatial 
order  of  accuracy  of  the  velocity  gradient  on  the  solution  is  examined. 

2.2  Force  Transformation:  Infinite  Plate  Spline 

The  structure  and  fluid  grids  are  linked  by  the  “infinite  plate  spline”  approach  im¬ 
plemented  by  Parker  (TIBI)  for  use  with  FLUENT.  This  implementation  has  already  shown 
accurate  results  for  several  transonic  flutter  cases  (1 161).  The  spline  determines  the  positions 
of  the  CFD  grid  points,  qa,  by  minimizing  the  bending  energy  of  a  thin  plate  passing  through 
the  structural  grid  points,  qs  (8).  The  CFD  grid  points  are  given  by 

{<?4  =  [S]  {gj  (17) 


where  the  spline  matrix,  S,  is 


5  =  Bas  [B-1  -  Bf'Rs  (  Rts  Bf 1  R.s  )  ^Rj Bf1]  1  +  Ra  (Rj BfxR,s)  1  RTS  B 


>T  E>  — 1 


(18) 


with  R  and  B  given  by 


Bn 


'  as 

Bs 

r2 


Xl 

y  i 

Zl 

X2 

2/2 

Z2 

Xn 

Vn 

Zn 

r2 
'  as 

ln( 

r2  ) 

'  as) 

{xai  “1“  ( Uai  Usj )  (^ai  %sj ) 

r2s  hi  (rf) 

(%si  Lsj )  T  {ysi  Usj)  T  (%si  ^sj ) 


(19) 


(20) 


The  subscript  a  in  Equations  C[8)  US  and  [20]  denotes  a  matrix  of  coordinates  from  the 
CFD  grid,  the  subscript  s  denotes  a  matrix  of  coordinates  from  the  structure  grid,  and 
the  subscript  as  denotes  a  matrix  containing  both  coordinates  (1 161).  The  forces  can  be 
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transformed  to  the  structure  grid  simply  by 


Fs  =  STFa  (21) 

One  drawback  of  using  FLUENT  is  that  the  forces  are  updated  at  every  time-step, 
but  not  every  sub-iteration  of  the  implicit  solver.  This  limitation  is  due  to  the  interface 
available  for  user  defined  functions  in  FLUENT.  This  is  unfortunate  since  the  forces  must 
be  time  lagged  for  use  in  the  structural  model.  This  has  a  detrimental  effect  on  the  temporal 
accuracy  that  can  be  achieved  with  an  implicitly  integrated  scheme,  or  the  commonly  used 
multi-stage  explicit  schemes.  In  addition,  the  fluid  solution  is  limited  to  first  order  temporal 
and  spatial  accuracy  for  solutions  with  deforming  meshes  in  FLUENT. 

2.3  Structural  Solver:  Modal  Analysis 

The  harmonic  response  of  a  structure  is  described  by 

(K  -  uj2M)  d  =  0  (22) 

where  K  is  the  stiffness  matrix,  M  is  the  mass  matrix,  a;  is  a  natural  frequency,  and  d  is  the 
vector  of  displacements  of  the  nodes.  Harmonic  vibration  is  a  state  of  balance  between  the 
forces  due  to  inertia  (u>2Md)  and  the  elastic  forces  (Kd)  (6).  Rewriting  equation  [221  makes 
it  easy  to  recognize  the  generalized  eigen-problem 


Kd  =  ui2Md  (23) 

with  the  eigenvalues  being  a  natural  frequency  squared  and  the  eigenvector  being  the  vector 
of  displacements,  or  modes,  which  satisfies  equation [231  for  the  corresponding  eigenvalue  (I)- 
This  research  uses  experimentally  determined  eigenvalues  and  eigenvectors  based  on  laser 
vibrometery  measurements  of  the  dynamic  response  of  the  MAV  in  air  and  in  vacuum. 

2.3.1  Structural  equation  of  motion  on  a  reduced  basis.  The  advantage  of  modal 
analysis  is  that  the  model  uses  a  few  of  the  low  frequency,  energy  containing  modes  as  a 
reduced-dimension  basis  for  the  structural  dynamics  problem.  The  system  can  be  solved 
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more  efficiently  since  the  orthogonal  basis,  <h,  allows  an  uncoupled  solution  to  the  structural 
system,  given  by 

<$>TK<5>d  =  w2<FTM4>d  (24) 

We  first  seek  to  diagonalize  the  stiffness  matrix,  K  so  that  the  structural  model  can 
be  efficiently  solved  as  a  decoupled  system.  Noting  that  K  is  a  real  symmetric  matrix,  and 
recalling  the  spectral  theorem,  which  states  that  a  real  symmetric  matrix  is  diagonalizable 
by  an  orthogonal  matrix  with  the  eigenvectors  of  the  matrix  we  wish  to  diagonalize  as  its 
columns (TI8]h  we  have  our  orthogonal  matrix,  4?  which  has  the  property: 

4*T<F  =  I  (25) 

which  gives  us  the  inverse  as  well,  4)-1  =  <f>T.  <J>,  sometimes  referred  to  as  the  modal  matrix, 
has  the  eigenvectors  of  the  stiffness  matrix  as  its  columns. 

Generally,  the  mode  shapes  are  mass  normalized  such  that 

drMd  =  1  (26) 

Also,  the  eigenvectors  are  mass-orthogonal  with  each  other  (shown  in  Appendix  [A),  that  is 

djMdj  =  0  (i^j)  (27) 

Usually  the  mass  matrix  is  a  known  quantity  that  is  developed  as  a  preliminary  step  in 
a  structural  analysis,  such  as  the  finite  element  method.  For  our  purposes  we  assume  a 
uniform  mass  distribution,  giving  a  mass  matrix  of  the  form 

M  =  {T^L)  I  (28) 

where  n  is  the  number  of  nodes,  and  mt0tai  is  the  total  mass  of  the  structure,  which  is  a 
lumped-mass  matrix  assuming  uniform  mass  distribution.  The  mass  normalization  of  the 
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eigenvectors  then  is  given  by 


dT  Md  =  mtotal  c2 


d 


norm 


C-^TTltotal  J 


d 


(29) 


In  Equation  [291  the  value  of  the  inner  product  before  normalization  is  represented  by  c2.  The 
measured  frequencies  of  the  MAV’s  harmonic  response  give  the  mass  to  stiffness  ratio,  but 
the  problem  of  determining  the  masses  of  the  structure’s  grid  points  is  under-constrained 
unless  n  modes  are  measured.  This  is  unlikely,  since  the  number  of  points  needed  for 
adequate  spatial  resolution  would  require  prohibitively  high  frequency  modes  to  be  excited. 
This  leaves  the  choice  of  weights  a  somewhat  open  question. 

Unlike  the  eigenvectors  from  a  mathematical  model,  the  experimentally  measured 
eigenvectors,  or  mode  shapes,  will  not  in  general  be  orthogonal  for  a  variety  of  reasons. 
The  structure  has  some  non-linearity  which  manifests  itself  as  coupling  between  the  mode 
shapes.  In  addition,  experimentally  exciting  a  single  mode  perfectly  is  difficult,  so  the 
measured  shape  will  have  some  small  components  of  the  other  modes  present. 

We  use  the  Gram-Schmidt  process  to  remove  the  non-orthogonal  components  from 
the  modes  (j!8j).  The  lowest  frequency  mode  is  retained  unaltered  (except  for  normaliza¬ 
tion),  and  the  parallel  components  of  each  subsequent  mode  are  removed.  The  algorithm  is 
straightforward  and  shown  in  Equations  [30]  and  [311  with  a  full  listing  in  Appendix  El 

For  all  of  the  j  vectors: 


1.  Subtract  off  the  parallel  components  of  previous  vectors 

3- 1 

Qj  =  aj  (qiao)  qi  (30) 

2—1 


2.  Normalize 


Qj 

*“13 


The  modal  matrix  is  now 


d*  = 


Ql  <?2  •  •  '  Qn 


(31) 


(32) 
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which  is  used  to  diagonalize  the  structure’s  stiffness  matrix 


[cu2]  = 

The  spectral  matrix  in  equation  134]  is  given  by 

0 

r  2i  0  w2  0 

M= 

•  • •  0  ul 


(33) 


(34) 


An  arbitrary  displacement  vector  d,  or  acceleration  vector  d  is  written  as  a  linear  combina¬ 
tion  of  the  modes. 

d  =  $z 
d  =  $z 

where  z  is  the  vector  of  modal  displacements.  Substituting  equation  [35]  into  the  undamped 
equation  of  motion  gives 

—  K &z  =  Fext  (36) 

Pre- multiplying  equation  [36]  by  4>T  gives 


<$>TM<S>z  -  <5>tK <S>z  =  <S>TFext 


(37) 


With  the  eigenvectors  mass-normalized,  substituting  the  spectral  stiffness  matrix  from  equa¬ 
tion  [33]  gives  the  modal  equation  of  motion  (with  no  damping) : 

z  +  [uj2]  z  =  $TFext  (38) 

Truncating  the  mode  set  to  the  first  m  modes  gives 

m 

d  ~  ^2  $iZi 
1=1 

m 

d  ~  ^2  $iZi 
i= 1 
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where  $  is  now  anmxn  matrix.  Equation  [38]  is  integrated  in  time  using  this  reduced  basis. 

Two  separate  formulations  for  the  time  integration  of  the  governing  equations  of  the 
structure  are  examined,  an  explicit  multi-stage  scheme,  and  a  semi-implicit  method. 

2.3.2  Multi-stage  Explicit.  We  use  a  four  stage  CFD  style  Runge-Kutte  scheme 
(source  listing  for  all  integration  routines  in  Appendix  [C)  to  integrate  the  structural  equa¬ 
tions  of  motion.  The  stage  coefficients  are 

1  1  1 

au  =  -,a  2  =  2>a3  =  -,«4  =  1  (39) 


The  second  order  system  given  by  equation  [38]  can  be  written  as  a  coupled  system  of 
first  order  ODE’s. 


d 

z 

z 

dt 

z 

$tF  -  u2z 

(40) 


which  can  be  integrated  by  an  arbitrary  multistage  scheme  by 


- 

k+ 1 

n 

- 

z 

— 

-  ak  (At) 

z 

z 

z 

4>tF  -  uj2z 

(41) 


where  k  denotes  the  stage  level  and  n  denotes  the  time  level,  the  final  stage  gives  the  n  +  1 
solution.  For  an  m  stage  scheme  at  least  second  order  accuracy  is  achieved  with  am~1  =  0.5 
and  am  =  1  (2). 


2.3.3  Newmark  -  semi-implicit.  The  Newmark  method  makes  substitutions  for 
displacements  and  velocities  at  the  n  +  1  time  level  with  Taylor  series  expansions  about  the 
n  time  level  allowing  the  acceleration  to  be  updated  using  only  n  time  level  values.  Then  a 
new  velocity  can  be  calculated  based  on  an  average  of  the  acceleration  between  the  n  and 
n  +  1  time  levels  (10). 

..„+i  _  K  ~  ^  (*n  +  (At)i"  +  (At)2F73) 

1  +  (Af)2u;2/6  1  j 


=  z"  +  ^  (zn  +  zn+1) 


•n+1  _  •  n 


(43) 
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zn+ 1  =zn  +  (A  t)zn  + 


(44) 


(2  zn  +  zn+1 ) 

The  performance  of  both  of  these  schemes  is  evaluated  in  section  13.41  by  integrating  the 
free-vibration  of  the  modes  in  response  to  a  displacement  perturbation,  and  then  they  are 
applied  to  the  full  aeroelastic  problem. 

2-4  Laser  Vibrometry 

Laser  vibrometry  consists  of  exciting  the  structure  and  then  measuring  the  phase  shift 
in  a  laser  beam  reflected  off  of  the  vibrating  structure  against  a  reference  beam  split  from 
the  same  source(JT7l).  The  unit  used  for  this  research  allows  three  measurement  heads  to 
be  positioned  to  determine  velocity  and  displacements  of  a  vibrating  structure  in  three 
dimensions.  However,  this  analysis  assumes  that  the  main  deformations  are  in  the  direction 
normal  to  the  chord  line  of  the  wing  so  only  one  measurement  head  is  needed. 

The  experimentally  measured  frequency  is  in  units  of  cycles  per  second  (Hz).  The  u ;  in 
equations [22] and [23] is  a  circular  frequency  (radians  per  second).  The  simple,  but  important 
conversion  is 

cu2  =  (27 r/)2  (45) 


The  air  around  the  vibrating  structure  adds  an  “apparent  mass”  effect  to  the  measured 
frequencies  (5]).  Because  the  flexible  wing  is  so  light  the  effect  of  the  apparent  mass  of  air  is 
significant  so  this  effect  needs  to  be  removed  from  the  measured  mode  frequencies  so  they 
can  be  used  in  the  aeroelastic  solution.  We  achieve  this  by  measuring  the  vibrations  of  the 
wing  in  a  vacuum  chamber  at  varying  ambient  pressure  levels.  This  data  is  used  to  develop 
a  frequency  correction  to  the  modes  used  in  the  model. 

The  vacuum  chamber  has  limited  volume  and  optical  access,  so  only  the  wing  of  the 
MAV,  separate  from  the  fuselage,  was  used  in  the  variable  ambient  pressure  tests.  This  is 
a  reasonable  way  to  develop  the  apparent  mass  correction  since  the  main  affect  will  be  on 
the  relatively  low  inertia  (mass)  membrane  wing  rather  than  the  significantly  more  massive 


16 


Figure  3:  Membrane  wing  on  shaker  table 


fuselage.  Changes  to  the  mode  shapes  due  to  the  wing  being  removed  from  the  fuselage  are 
not  critical  since  the  correction  is  made  to  the  frequency  of  the  mode,  not  its  shape. 

2-4-1  Wing.  Measurements  of  the  wing’s  motion  are  made  by  the  scanning  head 
at  specific  locations  over  the  wing’s  surface.  These  locations  become  the  “grid”  for  the 
structural  model.  The  damped  mode  shapes  are  measured  easily  in  air  by  exciting  the 
model  with  a  large  speaker.  Figure  [3]  shows  the  wing  attached  to  the  shaker  table  with  the 
speaker  used  to  excite  the  structure  in  the  background.  These  mode  shapes  are  shown  in 
Figure  31  Figure  51  shows  the  very  noisy  mode  shapes  at  the  higher  frequencies.  This  is 
likely  caused  by  removing  the  wing  from  the  fuselage  which  removes  the  constraint  on  the 
central  portion  of  the  wing’s  trailing  edge  where  it  attaches  to  the  rest  of  the  aircraft. 

The  undamped  modes  are  measured  in  a  vacuum  chamber  with  a  small  port  for  optical 
access,  and  the  model  is  excited  by  a  small  piezo-electric  patch  attached  to  the  wing.  Due  to 
the  restricted  optical  access  in  the  vacuum  chamber,  the  mode  frequencies  are  measured  by 
a  single  point  scan  on  the  wing,  chosen  in  a  location  that  is  likely  to  have  a  large  amplitude 
response  for  all  of  the  modes.  This  allows  frequency  information  to  be  measured,  but  no 
shape  information  can  be  obtained  from  a  single  measurement  point.  Therefore  the  mode 
shapes  themselves  are  measured  in  air  with  the  multi-point  scan  over  the  entire  wing. 

The  correction  to  the  frequency  of  the  damped  shape  (measured  in  air)  is  made  based 
on  the  measured  frequency  variation  of  the  vacuum  results  (3). 
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Figure  4:  Wing  only  measured  mode  shapes 
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Figure  5:  Membrane  wing  in  vacuum  chamber 
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Figure  6:  Effect  of  ambient  pressure  on  mode  frequencies 


The  apparent  mass  of  the  air  moving  around  the  vibrating  structure  causes  the  MAV’s 
wing  to  vibrate  at  a  lower  frequency  in  air  ([5]).  The  shift  is  most  easily  seen  on  the  power 
spectrum  in  Figure  |B]  for  the  second  mode.  The  spike  around  20  Hz  shifts  to  a  larger  ampli¬ 
tude  and  higher  frequency  as  vacuum  is  approached.  We  use  a  simple  linear  extrapolation 
(equation  [46])  to  remove  the  effect  of  air  on  the  wing’s  vibration.  Note:  the  validity  of  this 
adjustment  will  be  evaluated  in  a  separate  research  effort. 


f  adjust  —  f measured  (1  T  s') 


(46) 


where  s  is  based  on  the  frequency  shift  of  the  second  mode  shape,  shown  in  detail  in  Figure 
[7|  The  legend  of  Figure  0  also  shows  the  decrease  in  measured  damping  (2%  to  1.3%)  as 
the  air  is  removed  from  the  chamber. 


19.9  -  17.6 

s  =  - 

19.9 

The  higher  modes  are  not  excited  very  reliably  by  the  peizo-electric  patch  for  the  variable 
ambient  pressure  tests.  Figure  [6]  shows  the  lack  of  well  defined  amplitude  peaks  for  any  of 
the  higher  modes. 
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Figure  7:  Frequency  shift  of  2nd  mode  due  to  apparent  mass 


2.4.2  Full  Aircraft.  The  full  aircraft  is  shown  mounted  on  the  shaker  table  in 
Figure  SI  The  speaker  used  to  drive  the  structure  during  the  measurements  is  seen  in  the 
foreground  of  this  picture,  the  scanning  head  is  not  visible  in  this  view. 

The  measured  mode  shapes  are  shown  in  Figure  01  The  magnitude  of  the  displacements 
measured  by  the  laser  vibrometry  scanning  head  are  very  small,  but  the  magnitude  of  the 
mode  vector  is  arbitrary,  and  will  be  normalized  as  part  of  the  Gram-Schmidt  process. 

Table  [2]  shows  the  inner  product  of  each  eigenvector  with  every  other,  shown  in  equa¬ 
tion  5H 

di  ■  dj 

The  non-zero  off-diagonal  components  in  Table  0  show  that  the  raw  measurements  contain 
elements  of  the  other  mode  shapes.  This  is  understandable  since  we  are  exciting  a  real 
structure  with  an  imperfect  driver  (the  speaker),  perfectly  exciting  just  a  single  mode  is 
unlikely.  The  closer  the  frequencies  of  the  modes  are  to  each  other  the  more  likely  they  will 
be  coupled,  or  co-excited  by  the  speaker.  Notably  the  fourth  and  fifth  modes  are  very  close 
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Figure  8:  Full  aircraft  mounted  on  shaker  table 


Figure  9:  Full  aircraft  measured  mode  shapes 
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Table  2:  Inner  products  of  the  nine  measured  mode  shapes 


1.00 

0.37 

0.40 

-0.32 

-0.19 

-0.03 

0.04 

-0.11 

-0.02 

0.37 

1.00 

-0.11 

-0.37 

-0.11 

0.13 

0.08 

0.06 

-0.14 

0.40 

-0.11 
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-0.26 

-0.10 

0.07 

0.00 

0.09 

-0.06 

-0.32 

-0.37 

-0.26 

1.00 

0.83 

-0.09 

0.14 

-0.05 

-0.03 

-0.19 

-0.11 

-0.10 

0.83 

1.00 

0.12 

0.33 

0.07 

-0.12 

-0.03 

0.13 

0.07 
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together  in  frequency,  and  their  components  are  very  nearly  parallel  (inner  product  close 
to  unity).  The  non-orthogonal  components  of  the  measured  shapes  must  be  eliminated  so 
they  can  be  used  to  diagonalize  the  stiffness  matrix. 

All  of  the  off-diagonal  components  shown  in  Table  @  are  around  machine  zero  for  the 
modes  after  the  Gram-Schmidt  process.  The  measured  mode  shapes  are  shown  in  Figure  [9J 
and  the  result  of  applying  the  Gram-Schmidt  process  to  achieve  an  orthonormal  basis  are 
shown  in  Figure  10. 

In  the  ortho-normalization  process  nothing  is  subtracted  off  of  the  first  mode,  it  is  only 
normalized.  The  fourth  and  fifth  mode  are  very  close  together  in  frequency  and  so  there  is 
significant  overlap  between  the  two  measurements,  referring  back  to  Table  [21  we  see  that  the 
fourth  and  fifth  mode  are  nearly  coincident.  The  wing/fuselage  deformations  are  roughly 
the  same,  but  the  tail  orientation  is  slightly  different,  this  results  in  a  orthonormalized  shape 
that  is  significantly  different  from  the  measured  shape. 

The  power  spectrum  of  the  harmonic  response  of  the  MAV  is  shown  in  Figure  [XXI  The 
frequency  peak  of  the  second  mode  shape  is  highlighted  in  the  lower  portion  of  the  figure 
and  this  shape  is  illustrated  in  the  upper  portion. 

The  higher  mode  shapes  of  the  detached  wing  bear  little  resemblance  to  the  movement 
of  the  wing  on  the  full  aircraft.  This  is  expected  since  the  trailing  edge  is  completely 
unrestrained.  Referring  back  to  Figures  5J  and  (2  we  see  that  the  higher  modes  are  difficult 
to  excite,  and  do  not  have  very  well  defined  frequency  peaks.  This  is  the  reason  for  using 
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Figure  10:  Orthonormalized  eigenvectors 
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Figure  11:  Power  spectrum  of  harmonic  response,  and  second  mode  shape  visualization 
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Table  3:  Full  aircraft  frequency  adjustment 


Mode 

f ‘measured  (Hz) 

/ adjust  (Hz) 

1 

22.75 

25.38 

2 

25.63 

28.59 

3 

40.75 

45.46 

4 

47.38 

52.86 

5 

49.75 

55.50 

6 

68.50 

76.42 

7 

70.75 

78.93 

8 

89.13 

99.43 

9 

94.13 

105.01 

the  more  well  defined  frequency  shift  of  the  second  mode  to  make  the  adjustments  due  to 
ambient  pressure. 

Table 0] shows  the  adjusted  frequencies  for  the  full  aircraft.  This  frequency  adjustment 
uses  the  measurements  made  on  the  wing  only,  but  is  applied  to  the  frequencies  measured  for 
the  whole  aircraft  configuration.  This  is  reasonable  since  the  greatest  affect  of  the  apparent 
mass  is  felt  by  the  motion  of  the  light,  flexible  wing,  and  the  symmetric  bending  mode  is 
similar  for  both  configurations. 

2.5  Solution  Algorithm 

The  complete  solution  is  accomplished  following  the  procedure  illustrated  in  Figure  [121 
The  basic  aeroelastic  algorithm  consists  of  a  series  of  transformations.  Forces  calculated  by 
FLUENT  are  transformed  to  the  modal  model,  and  displacements  calculated  by  integrating 
the  modes  are  transformed  back  to  the  fluid  grid. 

All  of  the  transformation  matrices,  the  infinite  plate  splines,  and  the  modal  basis,  are 
calculated  and  stored  once  at  the  beginning  of  the  solution.  The  proper  multiplications 
are  then  made  at  each  time  step  by  calling  user  defined  functions  in  FLUENT.  In  this  way 
the  structure  gets  forces  from  the  fluid  through  the  infinite  plate  spline  at  every  time-step, 
and  the  fluid  grid  is  moved  according  to  the  structural  displacements,  also  transformed 
by  the  infinite  plate  spline  at  each  time-step.  Within  the  structural  model,  the  forces  are 
transformed  from  grid  points  to  modes,  which  are  integrated,  and  then  transformed  back 
to  grid  point  displacements. 
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Figure  12:  Solution  algorithm  flow  chart 


26 


III.  Results 


The  validation  of  the  aeroelastic  model  begins  by  comparing  the  computational  solu¬ 
tions  with  no  deformations  to  the  wind-tunnel  results  for  the  rigid  wing  MAV.  The 
static  deformations  of  the  MAV  are  then  incorporated  into  the  steady-state  fluid  solution, 
and  these  results  are  compared  to  the  time-averaged  force  measurements  made  in  the  wind 
tunnel  on  the  flexible  wing  MAV.  Finally,  the  dynamic  response  of  the  MAV  in  steady 
free-stream  and  gusty  conditions  is  investigated. 

3.1  Problem  domain  and  computational  grid 

The  MAV  is  modeled  inside  the  wind  tunnel  which  was  used  for  the  experimental 
investigation.  This  reduces  the  overall  size  of  the  grid  required  to  accurately  resolve  the  flow- 
field,  but  requires  the  volume  grid  to  be  re-triangulated  for  each  different  angle  of  attack. 
Figure  J31  shows  the  MAV  at  0°  a  in  the  wind  tunnel  test-section  used  for  this  research.  The 
test  section  dimensions  are  31”  x44”  x72”,  which  gives  a  span-to-tunnel  width  ratio  of  0.545 
(H).  The  comparison  to  experimental  results  is  also  more  meaningful  since  the  wall  effects 
are  not  removed  in  the  computational  solution  as  they  would  be  if  the  MAV  was  modeled  in 
free  flight.  The  blockage  correction  applied  in  (0)  amounts  to  less  than  0.1%  of  the  dynamic 
pressure,  while  the  cross  tunnel  variation  in  velocity  was  measured  to  be  on  the  order  of  2% 
of  the  free-stream  velocity,  so  this  correction  is  neglected  in  the  present  calculations. 

The  fluid  grid  has  approximately  640,000  tetrahedral  cells,  based  on  the  geometry 
definition  used  for  preliminary  CFD  investigation  presented  in  (115|.  The  cell  volumes  range 
from  approximately  1.8  x  10~12  to  4.1  x  10~4  m3.  Figure  HI  shows  the  surface  grid  of  the 
MAV,  and  figure  IS  shows  the  faces  of  the  interior  cells  in  a  slice  across  the  front  portion  of 
the  MAV. 

One  drawback  to  using  an  unstructured  grid  for  these  viscous  solutions  is  the  large 
number  of  points  required  on  the  MAV  surface.  The  unstructured  grid  spacing  along  the 
surface  must  be  on  the  order  of  the  boundary  layer  height  to  properly  resolve  this  near  wall 
feature  with  unit  aspect  ratio  pyramids.  A  structured  grid,  on  the  other  hand,  can  be  highly 
skewed  near  the  wall  to  achieve  much  a  higher  resolution  in  the  surface  normal  direction 
compared  to  the  surface  tangent  direction,  reducing  the  number  of  nodes  on  the  surface. 
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Figure  13:  MAY  in  wind  tunnel  test-section 
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Figure  14:  CFD  surface  grid  edges 


This  is  a  cost  for  any  analysis  using  unstructured  grids,  but  becomes  particularly  costly  in 
this  work  since  the  splines  connecting  the  movement  of  the  structure  and  CFD  grids  scale 
directly  with  the  number  of  surface  points.  For  these  grids  the  MAV  surface  has  on  the 
order  of  95k  faces  and  47k  points.  The  splines  connect  these  fluid  grid  components  to  the 
255  points  of  the  structure  “grid” . 

3.2  Steady  laminar  solutions 

This  series  of  solutions  is  an  alpha  sweep  of  steady,  non-deforming  solutions  at  10,  20 
and  30  mph  inlet  velocity.  The  free-streanr  properties  of  the  cases  are  shown  in  Table  15 
The  alpha  sweep  includes  -10,  -5,  -2,  2,  4-13,  15,  18,  and  20  degrees  angle  of  attack.  The 
three  tunnel  inlet  velocities  of  10,  20,  and  30  mph  were  chosen  to  match  the  nominal  inlet 
conditions  of  the  previous  wind  tunnel  investigation.  The  agreement  with  the  experimental 
data  is  slightly  improved  over  the  previous  investigation  in  (J15J)  by  increasing  the  resolution 
in  the  CFD  grid  near  the  wall,  especially  in  regions  of  high  curvature. 
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Figure  15:  CFD  volume  grid  slice 


Table  4:  Static  solution  ranges 


Angle  of  Attack 

Mach  Number 

Reynolds  Number 

Min 

-10 

1.3  x  10-2 

4.326  x  104 

Max 

20 

5.2  x  10~2 

1.730  x  105 

30 


Figure  16:  Lift  coefficient  variation  with  angle  of  attack  -  steady,  laminar 


The  lift  coefficient  for  steady-state  laminar  solutions  with  the  MAV  in  the  tunnel  is 
shown  in  Figure  [L6] .  These  solutions  were  run  with  all  rigid  surfaces  and  the  second  order 
upwind  momentum  discretization,  and  compare  well  to  the  experimental  data  for  the  rigid 
MAV  from  (7).  The  Reynolds  number  is  low  enough  that  a  completely  laminar  solution 
can  be  safely  assumed  up  until  significant  separation  occurs.  The  onset  of  stall  seems  to  be 
predicted  well  by  the  CFD  for  the  rigid  case. 

For  the  rigid  wing  MAV,  the  experimental  results  show  no  measurable  variation  with 
velocity  in  the  linear  region  of  the  lift  curve  (note  overlapping  error  bars  in  Figure  [16) .  The 
non-deforming  CFD  solution  also  shows  this  lack  of  variation  with  velocity,  and  predicts 
force  coefficients  closer  to  the  higher  velocity  wind  tunnel  data.  The  20  and  30  mph  cases 
run  in  the  wind  tunnel  show  the  expected  Reynolds  number  variation  in  stall  angle  of  attack: 
the  higher  Reynolds  number  case  separates  at  a  higher  angle  of  attack.  The  experiments 
measured  the  10  mph  case  to  have  a  delayed  stall  angle  of  attack  compared  to  the  two 
higher  Reynolds  number  cases,  which  is  contrary  to  the  behavior  generally  expected  of 
stalling  airfoils/aircraft.  As  mentioned  previously,  the  experimental  data  shown  in  Figure 
Bll  has  been  adjusted  slightly  for  wall  effects. 

As  the  MAV  stalls  with  increasing  angle  of  attack  the  measured  lift  and  the  computed 
lift  agree  closely  for  the  rigid  case,  up  to  massive  separation  around  13  degrees  angle  of 
attack,  at  which  point  the  laminar  CFD  solution  is  unsuited  to  modeling  the  detached, 
turbulent  eddies  and  free-shear  layers  off  the  suction  surface  of  the  vehicle. 

At  13  degrees  angle  of  attack  the  CFD  solutions  show  a  slight  bump  in  the  lift,  but  at 
this  point  the  flow  is  massively  separated  over  large  sections  of  the  wing.  Figure  (171  shows 
velocity  vectors  over  the  wing  on  a  plane  0.15  m  from  the  centerline  of  the  vehicle.  The 
laminar  solution  is  unreliable  for  this  type  of  flow-field  involving  large  turbulent  free-shear 
layers. 

Figure  [IS]  shows  the  drag  coefficient  for  the  non-deforming  CFD  solution  and  the 
rigid  wing  wind  tunnel  results.  The  computation  agrees  well  with  the  experimental  results 
in  the  lower  portion  of  the  “drag  bucket”  where  skin  friction  (wall  shear)  dominates  the 
drag.  This  is  another  indication  (aside  from  the  low  Reynolds  number)  that  the  solution  is 
predominantly  laminar,  since  the  wall  shear  is  much  greater  for  a  turbulent  solution. 
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Figure  17:  Velocity  vectors  at  50%  span  -  13  degrees  aoa,  20mph 

Though  there  is  negligible  difference  in  the  integrated  coefficients  in  the  linear  region 
for  the  different  velocities,  there  is  a  detectable  difference  in  the  flow  structure.  Figure  [19] 
shows  a  slight  laminar  separation  bubble  develops  around  60%  of  the  chord  for  the  9  degrees 
angle  of  attack,  10  mph  case,  which  is  not  present  in  the  higher  velocity  case  (Figure  [20]) . 
The  bubble  is  small  and  static  enough  that  the  steady  state  solution  still  converges  for  these 
low  velocity  runs.  The  effect  of  this  bubble  was  postulated  to  cause  the  waviness  in  the  lift 
coefficient  for  the  10  mph  wind  tunnel  data.  The  CFD  solutions  do  not  tend  to  support 
this,  given  the  size  of  the  bubble,  which  quickly  proceeds  to  full  separation  over  the  rear 
portion  of  the  wing.  The  more  likely  explanation  for  the  odd  low  speed  data  is  that  the 
forces  on  the  MAV  are  on  the  order  of  the  balance  resolution  for  this  low  velocity.  Some 
waviness  is  visible  in  the  stall  region  of  the  lift  curve  for  the  20  and  30  mph  cases,  but  these 
cases  do  not  show  the  separation  bubble  that  is  apparent  in  the  low  speed  case.  (pH). 
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Figure  19:  Velocity  vectors  at  50%  span  -  9  degrees  aoa,  lOrnph 
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Figure  20:  Velocity  vectors  at  50%  span  -  9  degrees  aoa,  20mph 


35 


Figure  21:  Lift  coefficient  -  1st  and  2nd  order  spatial  discretizations 


The  effect  of  the  spatial  order  of  accuracy  is  most  clearly  visible  in  the  failure  of  the 
first  order  scheme  to  predict  the  correct  stall  angle  of  attack.  Figure  [211  shows  the  continued 
linear  increase  in  lift  between  eight  and  ten  degrees  angle  of  attack  when  the  second  order 
solution  and  experimental  results  show  significant  stall  behavior. 

3.3  Statically  Deformed  Structure 

All  of  the  solutions  in  this  section  use  the  first  order  upwind  discretization  for  the 
momentum  equations  in  the  fluid  solver.  This  limits  the  applicability  of  the  model  to  the 
linear  region  of  the  lift  curve,  as  shown  previously  by  the  non-deforming  calculations. 

The  time  derivative  in  the  structural  model  is  neglected  for  these  solutions  to  cal¬ 
culate  the  static  deformations  in  steady  flight  conditions.  The  change  in  lift  due  to  the 
static  deformation  of  a  membrane  wing  results  from  a  balance  of  two  competing  types  of 
displacements.  The  trailing  edge  of  the  wing  can  move  up  higher  than  the  leading  edge, 
which  reduces  the  effective  angle  of  attack  (“washout”),  and  contributes  a  reduction  in  lift. 
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Figure  22:  Static  deformation  -  lOaoa,  20  mph 

A  membrane  wing  may  also  “billow”  in  the  wind,  which  changes  the  camber  of  the  wing, 
and  can  contribute  an  increase  in  lift  (120). 

The  deformed  shape  of  the  structure  at  10  degrees  angle  of  attack  is  shown  in  Figure 
[221  The  largest  amplitude  is  calculated  for  the  second  mode  shape  (symmetric  bending) 
over  the  full  angle  of  attack  range  tested  for  the  statically  deformed  MAV.  Since  the  loading 
at  this  flight  condition  is  symmetric,  this  is  the  mode  that  should  be  most  excited.  Figure  [231 
shows  the  amplitude  variation  with  angle  of  attack  for  the  20  mph  case.  One  troubling  result 
is  the  significant  amplitude  of  the  first  mode  shape  (anti-symmetric),  which  is  unexpected 
for  this  symmetric  loading  situation.  The  modal  amplitudes  increase  almost  linearly  with 
angle  of  attack,  as  would  be  expected  since  the  lift  curve  is  linear  in  this  region,  and  the 
modal  amplitudes  are  a  linear  function  of  the  force  (see  Figure  [23]).  The  modes  retain 
approximately  the  same  relative  magnitude  to  each  other  throughout  the  angle  of  attack 
range,  with  the  higher  frequency  modes  being  only  slightly  activated.  Even  the  largest 
deformations  are  only  on  the  order  of  10~5  m. 
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Figure  23:  Static  modal  amplitudes  -  20  mph 


The  significant  activation  of  the  anti-symmetric  mode  is  likely  the  result  of  the  or¬ 
thonormalization  performed  on  the  measured  mode  shapes.  The  modes  were  all  scaled  to 
unit  length  regardless  of  their  measured  amplitudes  and  a  uniform  mass  distribution  was 
assumed.  This  gives  greater  displacement  magnitudes  in  the  structure  due  to  the  loading 
than  would  be  present  in  the  actual  structure  with  a  more  physiclly  accurate  mass  matrix. 
The  loading  has  a  significant  anti-symmetric  component,  but  it  is  given  too  much  weight 
by  the  present  mode  normalization. 

A  more  accurate  mass  matrix  could  be  created  by  some  type  of  area  weighting  of  the 
mass  of  each  point,  or  a  reference  displacement  of  the  actual  MAV  could  be  used  to  scale 
the  modes  to  achieve  the  correct  displacement  magnitudes  in  the  computational  model. 

Figure  [23  shows  the  lift  calculated  using  a  statically  deformed  structure,  along  with 
the  wind  tunnel  measurements  for  the  flexible  wing  MAV.  The  computations  do  not  show 
the  velocity  dependence  in  the  lift  curve  that  the  wind  tunnel  experiments  measured.  Figure 
Ml  shows  that  the  lift  computed  for  20  and  30  mph  free-stream  velocities  are  nearly  identical. 
In  the  linear  region  of  the  lift  curve  (a  <  8°)  the  computed  solution  falls  within  the  error 
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Figure  24:  Lift  coefficient  -  statically  deformed  structure 


bars  of  the  30  mph  case,  but  fails  to  capture  the  stall  past  this  angle  of  attack.  This  failure 
was  shown  for  the  non- deforming  case  as  well. 

The  lift  computed  for  the  statically  deformed  structure  is  nearly  identical  to  that 
computed  for  the  rigid  structure.  Figure  [251  shows  the  nearly  indistinguishable  performance 
predicted  by  the  rigid  and  flexible  models.  This  behavior  is  somewhat  supported  by  the 
experimental  data.  Figure  1261  shows  the  overlapping  error  bars  of  the  rigid  and  flexible  wing 
for  the  10  and  20  mph  cases  in  the  linear  region  of  the  lift  curve.  The  30  mph  case  has  a 
measurable  difference  in  lift  even  in  this  region  that  the  computations  do  not  capture. 

3-4  Steady  Free-stream  -  Dynamic  Structure 

The  solutions  in  this  section  use  the  first  order  upwind  discretization  for  the  fluid 
momentum  equations  and  the  two  previously  discussed  structural  integration  schemes.  Pre¬ 
vious  studies  of  membrane  wings  have  shown  “self-initiated”  high  frequency  vibrations  in 
steady  flight  conditions.  This  type  of  behavior  is  evaluated  with  the  fully  dynamic  struc¬ 
tural  model.  The  time  step  size  for  the  dynamic  solution  must  be  chosen  based  on  the 
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Figure  27:  Runge-Kutta  -  O  (At)2 


highest  frequency  of  the  retained  modes,  which  in  this  case  is  the  ninth  mode  at  105  Hz. 
Figure  [27]  shows  the  response  of  this  mode  to  a  displacement  perturbation  with  a  time  step 
of  2.6559  x  10~3s,  or  about  4  steps  per  cycle.  The  numerical  error  with  a  time  step  this 
large  causes  artificial  damping  and  the  contribution  from  this  mode  is  nearly  non-existent 
after  only  50  iterations.  Even  with  ten  steps  per  cycle  (At  =  1.0624  x  10-3s)  the  numerical 
damping  of  the  highest  mode  is  slightly  noticeable.  With  twenty  time-steps  per  cycle  of  the 
highest  mode  and  higher  the  numerical  damping  is  negligible. 

Figure  [281  shows  the  response  of  the  ninth  mode  to  an  initial  displacement  perturbation 
with  no  external  forces  or  damping  integrated  by  the  Newmark  method.  Even  with  only 
four  time  steps  per  cycle  the  Newmark  method  is  able  to  integrate  the  mode  without  undue 
dissipation,  though  there  are  noticeable  variations  in  the  amplitude  at  this  resolution.  At 
higher  resolutions  the  two  schemes  behave  similarly. 

The  ability  of  either  scheme  to  produce  accurate  results  in  the  coupled  fluid-structure 
solution  is  more  restricted  than  the  simple  structure  integration  previously  presented  would 
indicate.  At  a  resolution  of  approximately  four  time  steps  per  cycle  of  the  highest  frequency 
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Figure  28:  Newmark  method  -  O  (At)2 


retained,  which  gives  around  twenty  samples  in  the  lowest  frequency,  the  explicit  scheme 
shows  severe  dissipation  (Figure  [27]),  but  the  Newmark  method  showed  some  promise  (Figure 
128)  at  capturing  the  response  in  the  unforced  structure  integration.  The  power  spectrum  of 
lift  is  shown  in  Figures  [291  and [301  for  a  time  step  of  0.002  seconds,  which  gives  approximately 
five  time  steps  per  cycle  of  the  highest  frequency  retained  in  the  structural  model.  Figure  [291 
shows  the  results  obtained  with  the  Runge-Kutta  scheme,  and  Figure  [30]  shows  the  results 
obtained  with  the  Newmark  method,  starting  from  a  statically  deformed  solution.  The 
vertical  lines  on  the  plots  are  the  frequencies  of  the  modes  in  the  structural  model.  Neither 
scheme  is  able  to  accurately  resolve  the  high  frequency  modes  as  in  Figures  [33]  and  [34] 

A  spurious  low-frequency  component  is  evident  in  both  cases,  most  likely  caused  by 
aliasing  of  fluid  physics  at  frequencies  higher  than  250  Hz  (since  the  time-step  is  below  the 
Nyquist  limit  of  the  structure’s  frequencies).  The  explicit  scheme  has  little  more  than  noise 
at  the  intermediate  to  high  frequencies,  while  the  implicit  scheme  produces  a  significant 
component  (at  about  95  Hz)  below  the  highest  frequency.  An  additional  spike  is  also  seen 
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Figure  30:  Cl  frequency  -  lOaoa,  At  =  0.002,  Newmark 


Hz 


Figure  31:  Cl  frequency  -  lOaoa,  At  =  0.001,  Runge-Kutta 

in  both  methods  between  the  second  and  third  mode,  the  Newmark  method  calculates  this 
component  to  be  larger  than  the  component  at  the  second  mode  shape. 

Reducing  the  time-step  to  0.001  s  does  not  noticeably  improve  the  performance  for 
either  scheme.  Figures  [311  and  [321  show  the  lift  variation,  which  is  completely  different  from 
the  resolved  solution  at  At  =  0.0005s  (shown  in  Figure  |33|).  The  aliasing  actually  seems  to 
be  worse  in  this  case  than  the  solution  with  At  =  0.002s.  As  mentioned  previously,  this  is 
most  likely  an  effect  of  the  non-linear  fluid  physics  being  under-sampled. 

The  variation  of  lift  coefficient  in  time  beginning  from  a  statically  converged  solution 
is  shown  in  Figure  [331  This  solution  is  integrated  with  the  explicit  scheme  examined  in 
section  12.3.21  using  a  time  step  of  5  x  10_4s,  which  is  approximately  twenty  time  steps  per 
cycle  of  the  highest  mode.  The  sharp  spikes  are  due  to  restart  errors  with  FLUENT,  but 
they  are  much  higher  frequency  than  the  relevant  variations  due  to  the  physics.  Figure  [331 
shows  the  power  spectrum  of  the  lift  coefficient  up  to  120  Hz.  The  vertical  lines  on  the 
graph  mark  the  natural  frequencies  of  the  modes.  At  this  time  step  we  see  significant  peaks 
at  the  frequencies  of  the  structural  model,  which  is  the  expeected  behavior.  This  shows  that 
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Figure  32:  Cl  frequency  -  lOaoa,  At  =  0.001,  Newmark 

the  characteristic  frequencies  of  the  fluid  behavior  are  significantly  higher  than  the  structure 
(since  the  structure  is  resolved  at  the  previously  presented  larger  time-steps) .  Surprisingly, 
significant  power  is  associated  with  the  higher  frequency  mode  around  100  Hz.  These  high 
frequency  fluctuations,  similar  to  those  reported  in  (111),  are  termed  “self-initiated”,  and  are 
due  to  small  velocity  and  pressure  perturbations  present  in  any  real  flow,  or  roundoff  errors 
present  in  a  computational  flow.  The  effect  of  large  flow  disturbances  is  further  evaluated 
in  Section  3.51 

The  frequency  adjustment  made  to  remove  the  affect  of  ambient  pressure  from  the 
measured  mode  shapes  shifted  the  frequencies  higher  than  those  measured  at  atmospheric 
pressure.  The  frequencies  of  the  response  would  be  expected  to  move  back  near  their 
measured  values  in  the  aeroelastic  solution,  but  they  actually  have  a  tendency  to  shift 
higher  than  the  frequencies  in  the  model.  This  is  especially  visible  in  the  higher  modes, 
shown  in  Figure  [Ml 

The  non-linear  coupling  between  the  modes  in  the  fluid-structure  interaction  requires 
that  all  of  the  modes  be  accurately  resolved.  Higher  resolution  than  would  be  expected 
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Figure  35:  10%  Gust  time  response  -  step  input,  lOaoa 


from  looking  simply  at  the  decoupled  structure  (Figures  BTT  and  121)  is  required  in  the  closely 
coupled  algorithm. 

3.5  Unsteady  Gust  Response 

The  aeroelastic  solution  method  is  used  to  evaluate  the  MAV’s  gust  response  and 
unsteady  performance  capabilities,  as  well  as  the  effectiveness  of  the  adaptive  washout 
mechanism  at  stabilizing  lift  production.  The  restart  errors  in  FLUENT,  which  were  only 
a  nuisance  in  the  previous  analysis,  are  now  more  troublesome,  since  the  solution  must  be 
restarted  after  applying  a  gust  at  the  inlet.  For  this  reason  the  first  two  time-steps  after 
the  gust  are  not  shown  in  the  subsequent  figures. 

Figure  [551  shows  the  response  of  the  flexible  and  rigid  MAV  in  a  free-stream  of  8.9408 
m/s  and  10  degrees  angle  of  attack  to  a  ten  percent  step  increase  in  the  free-stream  velocity. 

The  dynamic  structure  does  not  provide  any  noticeable  smoothing  of  the  lift  increase  over 
the  response  of  the  rigid  structure.  The  high  frequency  components  of  the  structure  are 
excited,  as  shown  in  Figure  [361  It  is  understandable  that  the  high  frequency  modes  of  the 
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Figure  36:  10%  Gust  frequency  response  (dynamic  structure)  -  step  input,  lOaoa 

MAV  would  play  an  important  role  in  the  vehicle’s  gust  response,  since  a  step  input  can  be 
thought  of  as  a  forcing  function  at  infinite  frequency. 

Since  the  MAV  will  experience  large  gusts  as  part  of  its  normal  operations  the  high 
frequency  vibrations  of  the  membrane  wing  are  significant  to  its  performance.  As  Figure 
[36]  shows,  the  high  frequency  modes  are  the  most  significant  part  of  the  lift  variation  in  the 
gust  response. 
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IV.  Conclusions 


The  aeroelastic  solution  method  presented  in  this  thesis  is  capable  of  modeling  the 
primary  behavior  of  the  membrane  wing  aerodynamics.  The  steady,  non-deforming 
laminar  solutions  predict  the  MAV’s  stall  angle  of  attack  when  a  2nd  order  upwind  scheme 
is  used,  but  fail  to  match  the  wind  tunnel  data  with  a  first  order  scheme.  The  dynamic 
performance  of  the  MAV  shows  significant  high  frequency  activity  even  under  steady  flight 
condition. 

This  flexible  wing  MAV  shows  mixed  results  in  achieving  all  of  the  claimed  advantages 
of  membrane  wings.  The  stall  performance  of  the  vehicle  can  not  be  properly  evaluated  with 
this  model,  since  the  first  order  discretization  fails  to  capture  the  stall  characterisitcs  of  the 
non-deforming  case.  The  “adaptive  washout”  mechanism  seems  to  provide  no  smoothing 
influence  on  the  lift  production  of  the  flexible  wing  over  the  rigid  wing  in  a  gust. 

4.1  Evaluation  of  the  Model 

The  modal  model  predicts  the  proper  linear  variations  in  modal  amplitudes  with 
increasing  MAV  angle  of  attack,  but  the  actual  size  of  the  wing  displacements  is  sensitive 
to  the  normalization  chosen  for  the  mode  shapes.  The  choice  of  normalization  has  a  direct 
scaling  effect  on  how  the  forces  translate  to  displacements  in  the  model,  and  is  shown  to  be 
an  area  in  need  of  future  improvement  and  validation. 

The  gust  response  is  difficult  to  evaluate  accurately  with  the  present  model  for  a 
variety  of  reasons.  The  scaling  of  the  modes  would  have  an  affect  on  the  time  response,  and 
the  poor  restart  capabilities  of  the  fluid  solver  make  the  first  two  iterations  after  a  restart 
unreliable. 

4-2  Design  Considerations 

The  performance  of  the  MAV  (both  computational  and  experimental)  indicates  that 
changes  to  the  wing  to  promote  “billowing”  type  deformations  rather  than  “washout”  type 
deformations  could  improve  the  lift  production  over  the  rigid  wing  in  the  linear  region  of 
the  lift  curve.  The  current  flexible  wing  design  slightly  under-performs  the  rigid  wing  at 
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nominal  angles  of  attack  (for  those  velocities  where  there  is  a  measurable  difference),  but 
shows  good  resistance  to  stall  over  the  rigid  wing  due  to  the  static  washout. 

The  unsteady  results  indicate  that  the  wing  is  too  stiff  to  provide  any  smoothing  in 
gusty  conditions,  the  dynamic  flexible  structure  response  being  nearly  as  step-like  as  that 
of  the  rigid  structure.  Other  MAV’s  that  have  been  reported  to  enjoy  this  smoothing  effect 
had  latex  rubber  as  the  inter-spar  material  rather  than  parachute  nylon,  which  might  have 
contributed  to  the  desirable  lift  smoothing. 

4-3  Considerations  for  Further  Analysis 

The  basic  method  presented  shows  promise,  but  there  are  significant  areas  for  improve¬ 
ment  in  this  aeroelastic  model.  Both  the  fluid  and  structural  components  have  significant 
shortcomings.  The  scope  of  the  method’s  applicability  can  also  be  extended  by  improving 
the  fluid  solver  to  handle  unsteady  detached  flows.  Also,  further  insight  into  the  effect  of 
variations  in  inter-spar  material  would  benefit  future  MAV  designers. 

There  are  two  approaches  to  improving  the  shortcomings  in  the  structural  model.  The 
first  would  be  to  keep  the  uniform  mass  matrix,  and  apply  some  type  of  scaling  to  the  mode 
amplitudes  based  on  the  measured  response  of  the  MAV.  The  second  is  keeping  the  present 
normalization  and  building  a  more  exact  mass  matrix  based  on  physical  measurements  or 
more  detailed  models  of  the  vehicle. 

In  the  first  approach  each  mode  could  be  normalized  by  a  ratio  of  its  amplitude  to 
the  largest  mode’s  amplitude.  A  scaling  by  the  ratio  of  u>2  of  the  mode  to  the  largest 
mode’s  lo2  could  be  used  to  correct  for  the  effect  of  larger  velocities  at  higher  frequencies 
giving  larger  measured  amplitudes  by  the  laser  vibrometry.  This  gives  an  effective  “modal 
mass”  (4).  The  problem  with  this  approach  is  the  unspecified  proportionality  of  a  unit 
displacement  of  the  largest  mode  to  a  the  magnitude  of  force  in  that  mode.  An  alternative 
would  be  constructing  a  more  exact  mass  matrix,  since  this  analysis  assumed  a  uniform 
mass  distribution.  The  assumption  is  not  greatly  detrimental  to  the  static  deformations  in 
this  case  since  the  wind  tunnel  data  showed  them  to  have  only  a  small  effect  on  the  lift. 
However,  the  effectively  increased  inertia  of  the  wing  (which  is  actually  lighter  than  the  rest 
of  the  aircraft)  will  cause  the  unsteady  response  of  the  structure  to  be  slower  due  to  the 
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increased  mass  to  stiffness  ratio.  An  improvement  might  involve  some  type  of  area  average 
along  with  a  weighting  based  on  what  parts  of  the  MAV  the  measurement  point  is  near 
(membrane,  spar,  leading  edge,  etc),  or  a  more  detailed  physical  model,  such  as  a  finite 
element  model.  The  utility  of  the  model  could  be  extended  most  directly  by  allowing  for 
sub-iteration  updates  (no  time  lagging)  of  the  force  transformation  between  the  fluid  and 
structure.  Extension  to  higher  order  spatial  accuracy  on  the  moving  grid  would  also  be 
advantageous. 

The  fluid  solver  needs  to  be  capable  of  higher  order  spatial  accuracy  to  capture  the  rel¬ 
evant  flow  physics  while  maintaining  a  reasonable  problem  size.  The  size  of  the  splines  used 
to  couple  the  fluid  and  structure  are  driven  by  the  number  of  points  on  the  vehicle’s  surface, 
and  can  quickly  make  doing  extensive  analysis  prohibitive  due  to  hardware  limitations. 

More  detailed  analysis  of  the  coupling  (if  any)  between  flow  separation,  shedding  and 
the  membrane  vibrations  can  be  performed  by  using  a  detached  eddy  simulation  (DES)  for 
the  flow  physics.  The  current  methods  being  applicable  only  to  mostly  attached  laminar 
flows. 

The  design  of  future  MAV’s  would  also  benefit  from  an  analysis  of  how  different  types 
of  inter-spar  material  affect  the  dynamic  response  of  the  wing,  most  importantly,  their  effect 
on  any  desirable  passive  control  properties  such  as  lift  smoothing  in  a  gust  response. 
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Appendix  A.  Mass-orthogonality  of  mode  shapes 

his  derivation  follows  that  presented  in  Appendix  C  of  (6).  The  harmonic  motion  of  a 
structure  with  no  viscous  damping  is  described  by 

(K  -  u?M )  d  =  0  (48) 


The  ith  and  jth  eigenpair  satisfy 


(. K  -  uffM)  dt  =  0 

and 

(K  -  uj]M )  dj  =  0 

respectively.  Premultiply  [49]  by  dj  and  by  dj,  giving 

djKdi  —  dJufMdi  =  0 

and 

djKdj  —  dfu'jMdj  =  0 
Since  both  K  and  M  are  symmetric 

djKdi  =  djKdj 
djMdi  =  djMdj 

and  subtracting  [52]  from  [5T]  gives 


For  ujf  ±  ufj 


(i 'jJj  —  cof)  djMdj  =  0 
dj  Mdi  =  0 


The  mode  shapes  of  different  frequencies  are  M  orthogonal. 


(49) 

(50) 

(51) 

(52) 

(53) 

(54) 
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Appendix  B.  Source  for  mode  shape  orthogonalization 

B.l  chkmalloc.c 

#include  <stdio.h> 

#include  <stdlib.h> 

#include  "chkmalloc .h" 

/*  this  routine  taken  from  C  programming  notes  by  Steve  Summit  */ 

void  *chkmalloc (size_t  sz) 

{ 

void  *ret  =  malloc(sz); 
if (ret  ==  NULL) 

{ 

fprintf (stderr ,  "Out  of  memory\n") ; 
exit (EXIT_FAILURE) ; 

> 

return  ret ; 

> 

B.2  getwords.c 

#include  "getwords.h" 

/*  this  routine  taken  from  C  programming  notes  by  Steve  Summit  */ 

int  getwords (char  *line,  char  *words  [] ,  int  maxwords) 

{ 

char  *p  =  line; 
int  nwords  =  0; 

while (1) 

{ 

while (isspace(*p) ) 

p++; 


if (*p  ==  ’\0’) 
return  nwords ; 

words  [nwords++]  =  p; 

while (! isspace (*p)  &&  *p  !=  ’\0’) 

p++; 


if (*p  ==  ’\0’) 
return  nwords ; 

*p++  =  ’\0’ ; 
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if (nwords  >=  maxwords) 
return  nwords ; 

> 

> 

B.3  ortho-norm. c 

/*  program  to  create  an  orthonormal  basis  of  eigenvectors  from 
measured  mode  shapes  */ 

#include  " ortho-norm. h" 

#include  "getwords.h" 

#def ine  MAXLINE  200 
#def ine  MAXWORDS  20 

int  main(int  argc,  char  *argv[]) 

{/*  first  argument  is  a  text  file  with  n,  and  m,  and  then  the  xy 
coordinates  of  the  nodes,  the  second  argument  is  the  text  file 
with  the  modes  as  its  columns  */ 

RKIND  **A;  /*  storage  for  the  modes  */ 

RKIND  **B;  /*  storage  for  the  amount  subtracted  off  in  Gramm-Schmidt  */ 

RKIND  *y;  /*  y  coordinate  of  nodes  */ 

RKIND  *x;  /*  x  coordinate  of  nodes  */ 

RKIND  *scale;  /*  for  normalizing  each  eigenvector  */ 

RKIND  *omega;  /*  natural  frequency  */ 

RKIND  **ortho;  /*  check  for  mutual  orthogonality  of  eigenvectors  */ 

IKIND  *n;  /*  rows  */ 

IKIND  *m;  /*  columns  */ 

RKIND  *mass;  /*  total  mass  of  the  structure  */ 
int  i ; 
int  j; 
int  k; 

FILE  *A_file;  /*  file  with  measured  shapes  as  columns  */ 

FILE  *valu_file;  /*  file  with  some  values  in  it  */ 

FILE  *ortho_tab;  /*  output  the  results  */ 

FILE  *w_file;/*  file  with  the  measured  freqs  in  Hz  */ 

FILE  *w_file_out; 
char  line [MAXLINE] ; 
int  ac ; 

char  *av [MAXWORDS]  ; 

n  =  (IKIND  *) chkmalloc (  sizeof (IKIND)  ); 
m  =  (IKIND  *) chkmalloc (  sizeof (IKIND)  ); 
mass  =  (RKIND  *) chkmalloc (  sizeof (RKIND)  ); 
scale  =  (RKIND  *) chkmalloc (  sizeof (RKIND)  ); 

if  (argc  !  =4){printf  ("Usage :  \n  7»s  value_f  ilename  A_filename  w_f  ile\n" ,  argv  [0] ) ;  exit(l' 


56 


else{ 

A_file  =  f open(argv  [2] , "r") ; 

if (A_f ile==NULL) {printf ("Can’t  open  °/0s  for  read.\n",  argv[2]);  exit(l);} 
valu_file  =  f open(argv [1] , "r") ; 

if (valu_file==NULL){printf ("Can’t  open  °/0s  for  read.\n",  argv[l]);  exit(l);} 
w_file  =  f open(argv  [3] , "r") ; 

if (w_f ile==NULL) {printf ("Can’t  open  °/0s  for  read.\n",  argv[3]);  exit(l);} 
w_file_out  =  fopenCstiffness.dat",  "w"); 

> 

/*  read  in  the  number  of  modes,  the  number  of  nodes,  and  the  coordinates  */ 

fgets(line,  MAXLINE,  valu_file) ; 

ac  =  getwords (line ,  av,  MAXWORDS-1) ; 

n[0]  =  atof(av[0]); 

m[0]  =  atof(av[l]); 

mass  [0]  =  atof(av[2]); 

y  =  (RKIND  *) chkmalloc (  m [0] *sizeof (RKIND)  ); 
x  =  (RKIND  *)chkmalloc(  m [0] *sizeof (RKIND)  ); 
omega  =  (RKIND  *) chkmalloc (  n  [0] *sizeof (RKIND)  ); 
i  =  0; 

while (f gets (line,  MAXLINE,  w_file)  !=  NULL){ 
ac  =  getwords (line ,  av,  MAXWORDS-1); 
omega [i]  =  atof (  av[0]  ); 

omega [i]  =  TW0*M_PI*omega[i] *TW0*M_PI*omega[i] ; 
fprintf  (w_f  ile_out,  "°/0l .  12e\n" ,  omega  [i] )  ; 
i++; 

> 

i  =  0; 

while (f gets (line,  MAXLINE,  valu_file)  !=  NULL){ 
if  (i<m[0] ){ 

ac  =  getwords (line ,  av,  MAXWORDS-1); 
x[i]  =  atof(av[0]); 
y  [i]  =  atof  (av[l]) ; 

}else{ 

printf  ("Wrong  number  of  nodes.  Specified:  °/0i,  Reached:  %i\n"  ,m[0]  , i) ; 
exit  (1) ; 

} 

i++; 

> 

A  =  (RKIND  **) chkmalloc (  n[0] *sizeof (RKIND)  ); 

B  =  (RKIND  **) chkmalloc (  n[0] *sizeof (RKIND)  ); 
ortho  =  (RKIND  **) chkmalloc (  n  [0] *sizeof (RKIND)  ); 
for(i=0;  i<n[0] ;  i++){ 

/*  A  and  B  are  (n  x  m) ,  ortho  is  only  (n  x  n)  */ 

A [i]  =  (RKIND  *) chkmalloc (  m [0] *sizeof (RKIND)  ); 

B  [i]  =  (RKIND  *) chkmalloc (  m  [0] *sizeof (RKIND)  ); 
ortho  [i]  =  (RKIND  *) chkmalloc (  n [0] *sizeof (RKIND)  ); 
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/*  the  mode  file  for  aeroelastic . c  has  the  eigenvectors  as  its  columns  */ 

j  =  0; 

while (f gets  (line,  MAXLINE,  A_file)  !=  NULLH 
ac  =  getwords (line ,  av,  MAXW0RDS-1) ; 
if  (  ac  !=  n[0]  ){ 

printf  ("Wrong  number  of  modes  found  in  °/0s .  \n" , argv  [2] ) ; 
exit (1) ; 

}else-[ 

for(i=0;  i<n[0] ;  i++){ 

A  [i]  [ j ]  =  (RKIND) atof (av  [i] ) ; 

> 

} 

j  =  j  +  i; 

> 

if (  j  !=  m[0]  ){ 

printf  ("Wrong  number  of  points  found  in  °/0s  .  \n" ,  argv  [2] ) ; 
exit (1) ; 

> 

f close (A_f ile) ; 

/*  print  out  the  modes  for  plotting  */ 

A_file  =  fopenCmode_exp.dat",  "w"); 
for (j=0;  j <m [0] ;  j++){ 

f  printf  (A_f  ile,  "%g  %g  ",x[j],  y[j]); 
for(i=0;  i<n[0] ;  i++){ 

fprintf  (A_f  ile,  "°/,g  " ,  A [i]  [ j ] ) ; 

} 

fprintf (A_f ile ,  "\n"); 

> 

f close (A_f ile) ; 

/*  normalize  each  row,  first  calculate  the  inner  product  of  the 
vector  with  itself,  then  divide  all  of  the  entries  by  the 
square  root  of  the  result  */ 
for(i=0;  i<n[0] ;  i++){ 

scale  [0]  =  inner  (A  [i]  ,  A[i],  m)  ; 
scale [0]  =  sqrt (scale  [0] ) ; 
f or (j=0;  j <m [0] ;  j++){ 

A  [i]  [j]  =  A[i]  [j]  /scale  [0]  ; 

} 

> 

/*  Check  and  see  if  the  measured  modes  are  close  to  orthogonal  with 
each  other,  then  print  it  out  for  a  latex  table. 

The  inner  product  of  the  i~th  vector  with  the  j~th  vector,  j  != 
i,  should  be  close  to  zero  */ 


58 


ortho_tab  =  fopen("ortho_exp_tab.tex" , "w") ; 
if(ortho_tab  ==  NULL){ 

printf ("Couldn't  open  ortho_exp_tab.tex  for  write. \n"); 
exit (1) ; 

> 

f printf  (ortho_tab,  "Wbegin-ttabularM  |  c  |  c  |  c  |  c  |  c  I  c  |  c  |  c  |  c  |  c}  \\hline  \n"); 
for(i=0;  i<n[0] ;  i++){ 
f or (j=0;  j <n [0] ;  j++){ 

ortho  [i][j]  =  inner (A [i] , A [j] ,m) ; 
if (j  ==  n [0] -1){ 

if  (  ortho[i][j]  <  0)  f  printf  (ortho_tab ,  "  °/01.2f  ",  ortho  [i]  [j] ) ; 
else  f printf (ortho_tab ,  "  %1.2f  ",  ortho [i] [j] ) ; 

}else-[ 

if  (  ortho[i][j]  <  0)  f  printf  (ortho_tab ,  "  °/01.2f  &  ",  ortho  [i]  [j] ) ; 
else  f printf (ortho_tab ,  "  %1.2f  &  ",  ortho  [i]  [j] ) ; 

> 

} 

f printf (ortho_tab ,  "\\\\  \\hline  \n") ; 

> 

f printf  (ortho_tab,  "\\end{tabular}\n") ; 
f close (ortho_tab) ; 

/*  Of  course  they  aren’t  orthogonal,  so  use  Gramm-Schmidt  to 
subtract  off  the  non-orthogonal  parts  */ 
for(i=l;  i<n[0] ;  i++){/*  loop  over  the  modes  */ 

for(k=i;  k>0;  k — ){/*  loop  over  the  previously  orthogonalized  vectors  */ 
scale  [0]  =  inner(A[k-l] ,A[i] ,m) ; 
for(j=0;  j<m[0] ;  j++){ 

B[i][j]  =  B[i][j]  +  scale  [0]  *A  [k—  1]  [j]  ;  /*  store  the  amount  we 
are  subtracting 
off,  to  look  at 
latter  */ 

A  [i]  [j]  =  A  [i]  [j]  -  scale  [0]  *A  [k-1]  [j]  ; 

> 

} 

/*  now  normalize  the  new  vector  */ 
scale [0]  =  inner (A  [i] ,  A[i] ,  m) ; 
scale  [0]  =  sqrt (scale  [0] ) ; 
f or (j=0;  j <m [0] ;  j++){ 

A  [i]  [j]  =  A[i]  [j] /scale  [0]  ; 

} 

> 

/*  print  it  out  just  to  see  if  we  did  it  right,  should  be  the  identity  */ 
ortho_tab  =  f open("ortho_GS_tab .tex" , "w") ; 
if(ortho_tab  ==  NULL)-( 

printf ("Couldn’t  open  ortho_GS_tab.tex  for  write. \n") ; 
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exit  (1) ; 


> 

fprintf  (ortho_tab,  "Wbegin-ttabularM  |  c  |  c  |  c  |  c  |  c  |  c  |  c  |  c  |  c  |  c}  \\hline  \n") 
for(i=0;  i<n[0] ;  i++){ 
f or (j=0;  j <n [0] ;  j++){ 

ortho  [i][j]  =  inner (A [i] , A  [j] ,m) ; 
if (j  ==  n  [0] -1){ 

if  (  ortho[i][j]  <  0)  fprintf  (ortho_tab ,  "  °/01.2f  ",  ortho  [i]  [j] ) ; 
else  fprintf (ortho_tab,  "  %1.2f  ",  ortho [i] [j] ) ; 

}else-[ 

if  (  ortho[i][j]  <  0)  fprintf  (ortho_tab ,  "  °/01.2f  &  ",  ortho  [i]  [j] )  ; 
else  fprintf (ortho_tab ,  "  %1.2f  &  ",  ortho  [i]  [j] ) ; 

> 

} 

fprintf (ortho_tab ,  "\\\\  \\hline\n") ; 

> 

fprintf  (ortho_tab,  "\\end-[tabular}-\n") ; 
f close (ortho_tab) ; 

/*  print  out  the  modes  for  plotting  */ 

A_file  =  fopenCmode_ortho.dat",  "w")  ; 
for (j=0;  j <m [0] ;  j++){ 

fprintf  (A_file,  "%g  °/.g  ",x[j],  y[j]); 
for(i=0;  i<n[0] ;  i++){ 

fprintf  (A_file,  "°/og  " , A [i]  [ j ] ) ; 

} 

fprintf (A_file ,  "\n"); 

> 

f close(A_f ile) ; 

/*  print  out  the  subtracted  components  for  plotting  */ 

A_file  =  fopenCmode_diff.dat",  "w"); 
f or (j=0;  j <m [0]  ;  j++){ 

fprintf  (A_f  ile,  "%g  °/,g  ",x[j],  y[j]); 
for(i=0;  i<n[0] ;  i++){ 

fprintf  (A_f  ile,  "°/,g  " , B [i]  [ j ] ) ; 

} 

fprintf (A_f ile ,  "\n"); 

> 

f close (A_f ile) ; 

/*  output  the  file  needed  for  the  structural  model  */ 
printf("Mass  norm:  %g,  °/od,  °/0g\n" , mass  [0]  ,m  [0]  ,  sqrt  (mass  [0] /m  [0] )) ; 
A_file  =  fopenCmodes.dat",  "w"); 
for (j=0;  j  <  m[0]  ;  j++){ 
for(i=0;  i  <  n[0] ;  i++){ 

fprintf  (A_f  ile ,  "°/01.12e  ",  A  [i]  [j] /sqrt  (mass  [0] /m  [0] )  ); 
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} 

fprintf (A_f ile, "\n") ; 

> 

f close (A_f ile) ; 
return(O) ; 

> 

RKIND  inner (RKIND  *vl,  RKIND  *v2,  IKIND  *len) 

{ 

RKIND  p  =  0; 
int  i ; 

for(i=0;  i<len[0] ;  i++){ 
p  =  p  +  vl  [i]  *v2  [i]  ; 

> 

return (p) ; 
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Appendix  C.  Source  for  structural  model  time  integration 


The  structural  model  is  integrated  in  time  by  an  explicit  multistage  Runge-Kutte  scheme 
as  well  as  a  second  order  semi-implicit  scheme.  A  sample  input  deck  for  the  program 
is  shown  below,  hash  marks  denote  comments. 


#  mode . input 

# 

#  An  input  deck  for  time  integration  of  the  modal  equations  of  motion 

#  generally  keyword  followed  by  value(s),  can  be  in  any  order,  but  all  of 

#  them  are  usually  required 

########### 

##  Files  ## 

########### 

#  the  mode  file  contains  the  orthonormal  mode  shapes  in  columns 
mode-file  modes.dat 

#  the  2D  location  of  each  of  the  sample  points 
xy-file  xy.dat 

#  this  file  contains  the  natural  frequencies  squared  in  rad/s,  the  first  value 

#  is  the  frequency  associated  with  the  shape  in  the  first  column,  etc 
spectral  stiffness.dat 

################ 

##  Parameters  ## 

################ 

#  number  of  modes 
n  9 

#  number  of  points 
m  255 

#  integration  method 
method  rk 


#  step  size 

#  four  steps  per  cycle  of  the  highest  freq 
#h  2 . 6559014129395500e-3 

#  ten  steps  per  cycle  of  the  highest  freq 
#h  1 . 06236056517582e-3 

#  fifteen  steps  per  cycle  of  the  highest  freq 
#h  7 . 08240376783881e-4 

#  twenty  steps  per  cylce  of  the  highest  freq 
h  5 . 31 1802825879 le-4 

#  total  number  of  time  steps  to  take 
iterations  2000 

#  write  data  this  many  iterations 
report  1 
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#  there  should  be  n  values  on  the  line  after  perturb,  the  integer  after 

#  perturb  is  0  for  displacement,  1  for  velocity 
perturb  0 

-1.0  0.8889  -0.7778  0.6667  -0.5556  0.4444  -0.3333  0.2222  -0.1111 

C.l  chkmalloc.c 

#include  <stdio.h> 

#include  <stdlib.h> 

#include  "chkmalloc .h" 

/*  this  routine  taken  from  C  programming  notes  by  Steve  Summit  */ 

void  *chkmalloc (size_t  sz) 

{ 

void  *ret  =  malloc(sz); 
if (ret  ==  NULL) 

{ 

fprintf (stderr ,  "Out  of  memory\n") ; 
exit (EXIT_FAILURE) ; 

> 

return  ret ; 

> 

C.2  getwords.c 

#include  "getwords.h" 

/*  this  routine  taken  from  C  programming  notes  by  Steve  Summit  */ 

int  getwords (char  *line,  char  *words  [] ,  int  maxwords) 

{ 

char  *p  =  line; 
int  nwords  =  0; 

while (1) 

{ 

while (isspace(*p) ) 

p++; 


if (*p  ==  ’\0’) 
return  nwords ; 

words  [nwords++]  =  p; 

while ( ! isspace (*p)  &&  *p  !=  ’\0’) 

p++; 
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if (*p  ==  ’\0’) 
return  nwords ; 

*p++  =  ’\0’ ; 

if (nwords  >=  maxwords) 
return  nwords ; 

> 

} 

C.  3  input,  c 

#include  " input. h" 

#def ine  MAXLINE  200 
#def ine  MAXWORDS  20 

void  alloc_input (INPUT  *in) 

{ 

in->n  =  (IKIND  *) chkmalloc (sizeof (IKIND) ) ; 
in->m  =  (IKIND  *) chkmalloc (sizeof (IKIND) ) ; 
in->its  =  (IKIND  *) chkmalloc (sizeof (IKIND) ) ; 
in->h  =  (RKIND  *) chkmalloc (sizeof (RKIND) ) ; 
in->report  =  (IKIND  *) chkmalloc (sizeof (IKIND) ) ; 
in->method  =  (IKIND  *) chkmalloc (sizeof (IKIND) ) ; 
in->n[0]  =  0; 
in->m[0]  =  0; 
in->its[0]  =  0; 
in->h[0]  =  0; 


void  get_input (FILE  *deck,  INPUT  *in,  RKIND  ***A,  RKIND  ***x,  \ 
RKIND  ***z,  RKIND  ***res,  RKIND  **omega) 

{ 

char  line [MAXLINE] ; 
int  ac ; 

char  *av [MAXWORDS]  ; 

IKIND  i , j ; 
alloc_input (in) ; 
j=0; 

while (  f gets (line,  MAXLINE,  deck)  !=  NULL  H 
if(*line  ==  ’#’) 
continue ; 

ac  =  getwords (line ,  av,  MAXW0RDS-1) ; 
if  (ac  ==  0) 
continue ; 
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if (  strcmp(av[0] ,  "mode-file")  ==  0  ){ 
in->mode_f ile  =  fopen(  av[l],"r"  ); 
if (in->mode_f ile  ==  NULL)  file_error(  av[l]  ); 

} 

if (  strcmp(av[0] ,  "xy-file")  ==  0){ 
in->xy_file  =  fopen(  av[l] ,  "r"); 
if (in->xy_f ile  ==  NULL)  file_error(  av[l]  ); 

} 

if (  strcmp(  av[0] ,  "spectral")  ==  0){ 
in->spectral  =  fopen(  av[l] ,  "r"); 
if (in->spectral  ==  NULL)  file_error(  av[l]  ); 

} 

if (  strcmp(  av[0]  ,  "n")  ==  0){ 
in->n[0]  =  atoi(  av[l]  ); 

} 

if (  strcmp(  av[0] ,  "m")  ==  0){ 
in->m[0]  =  atoi(  av[l]  ); 

} 

if (  strcmp(  av[0]  ,  "h")  ==  0){ 
in->h[0]  =  atof(  av[l]  ); 

} 

if (  strcmp(  av[0] ,  "iterations")  ==  0){ 
in->its  [0]  =  atof (  av[l]  ); 

} 

if (  strcmp(  av[0] ,  "report")  ==  0){ 
in->report [0]  =  atoi(  av[l]); 

} 

if (  strcmp(  av[0]  ,  "method")  ==  0){ 
if (  strcmp(  av[l] ,  "rk"  )  ==  0  ){ 
in->method [0]  =  0; 

> 

if(strcmp(  av[l] ,  "semi-implicit")  ==  0){ 
in->method  [0]  =  1; 

> 

if(strcmp(  av[l] ,  "third"  )  ==  0  ){ 
in->method [0]  =  2; 

> 

if(strcmp(  av[l] ,  "fourth"  )  ==  0  ){ 
in->method [0]  =  3; 

> 

if(strcmp(  av[l] ,  "fifth"  )  ==  0  ){ 
in->method [0]  =  4; 

> 

if(strcmp(  av[l] ,  "sixth"  )  ==  0  ){ 
in->method [0]  =  5; 

> 
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if(strcmp(  av[l] ,  "eighth"  )  ==  0  ){ 
in->method [0]  =  6; 

> 

if(strcmp(  av[l]  ,  "coupled-implicit-2"  )  ==  0  ){ 
in->method [0]  =  7; 

> 

if(strcmp(  av[l]  ,  "coupled-implicit-3"  )  ==  0  ){ 
in->method [0]  =  8; 

> 

if(strcmp(  av[l] ,  "coupled-implicit-4"  )  ==  0  ){ 
in->method [0]  =  9; 

> 

if(strcmp(  av[l] ,  "newmark"  )  ==  0  ){ 
in->method [0]  =  10; 

> 

if(strcmp(  av[l] ,  "var-accel-2"  )  ==  0  ){ 
in->method [0]  =  11; 

> 

if(strcmp(  av[l] ,  "rk-3"  )  ==  0  ){ 
in->method [0]  =  12; 

> 

} 

if (  strcmp(  av[0] ,  "perturb")  ==  0){ 
if (in->n  [0]  ==  0){ 

printf ("Declare  number  of  modes  before  perturbations\n") ; 
exit (1) ; 

> 

f gets (line,  MAXLINE,  deck); 

ac  =  getwords (line ,  av,  MAXW0RDS-1) ; 

if (  ac  !=  in->n[0]  ){ 

printf ("Wrong  number  of  perturbation  values  specif ied\n") ; 
exit (1) ; 

> 

in->perturb  =  (RKIND  *) chkmalloc(in->n [0] *sizeof (RKIND) ) ; 
for(i=0;i<in->n[0] ;i++)  in->perturb [i]  =  atof(  av[i]  ); 

} 

j++; 

> 

/*  done  parsing  the  input  deck  */ 

if (  (in->mode_f ile  ==  NULL)  ||  (in->xy_file  ==  NULL)  ){ 

printf ("Must  specify  mode  file  and  grid  file  in  input  deck\n") 
exit (1) ; 

> 

(*A)  =  (RKIND  **) chkmalloc (  in->n [0]  *sizeof (RKIND)  ); 
for(i=0;  i  <  in->n[0] ;  i++) 

(*A) [i]  =  (RKIND  *)chkmalloc(  in->m [0] *sizeof (RKIND)  ); 
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(*x)  =  (RKIND  **) chkmalloc (  DIM*sizeof (RKIND)  ); 
for(i=0;  i  <  DIM;  i++) 

Ox)  [i]  =  (RKIND  *) chkmalloc(  in->m [0] *sizeof (RKIND)  ); 

Oomega)  =  (RKIND  *) chkmalloc (  in->n [0] *sizeof (RKIND)  ); 

( *z)  =  (RKIND  **) chkmalloc (  HIST*sizeof (RKIND)  ); 

Ores)  =  (RKIND  **) chkmalloc (  3*sizeof (RKIND)  ); 
f or (i=0;  i  <  HIST;  i++){ 

(*z)  [i]  =  (RKIND  *)chkmalloc(  in->n [0] *sizeof (RKIND)  ); 

Ores)  [i]  =  (RKIND  *) chkmalloc (  in->n [0] *sizeof (RKIND)  ); 

> 

/*  read  the  grid  coordinates  */ 
i  =  0; 

while (  f gets (line,  MAXLINE,  in->xy_f ile)  !=  NULL  ){ 
ac  =  getwords (line ,  av,  MAXW0RDS-1) ; 
if (i<in->m[0] ){ 

Ox)  [0]  [i]  =  atof  (  av[0]  ); 

Ox)  [1]  [i]  =  atof  (  av  [1]  ); 

}else-[ 

printf("Too  many  points  in  grid  f ile\n") ; 
exit  (1) ; 

} 

i++; 

> 

j  =  0; 

while (  f gets (line,  MAXLINE,  in->mode_f ile)  !=  NULLM 
ac  =  getwords (line ,  av,  MAXW0RDS-1) ; 
if (ac  ==  in->n[0]M 

f or (i=0 ; i<in->n  [0] ;i++M 
(*A) [i]  [j]  =  atof (  av [i]  ); 

> 

}else-(  printf ("Wrong  number  of  points  on  the  line\n");  exit(l);  } 

j++; 

> 

j  =  0; 

while (  f gets (line,  MAXLINE,  in->spectral)  !=  NULLM 
ac  =  getwords (line ,  av,  MAXW0RDS-1) ; 
if (j<  in->n [0]  ) 

Oomega)  [j]  =  atof  (  av[0]  ); 
else{ 

printf ("Wrong  number  of  natural  frequencies  given\n"); 
exit  (1) ; 

} 

j++; 

> 

> 
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void  f ile_error (char  ^filename) 

{ 

printf  ("Error  with  10  attempt  to  °/0s\n" ,  filename); 
exit  (1) ; 


C-4  modal-  eom.c 

#include  "modal-eom.h" 

#def ine  MAXLINE  200 
#def ine  MAXWORDS  20 

int  main(int  argc,  char  *argv[]) 

{/*  orthonormalize  the  modes  before  using  them  in  this  program  */ 

FILE  *in_file; 

FILE  ^response; 

FILE  *mode_out; 

IKIND  j ,r; 

RKIND  **A; 

RKIND  **x;  /*  x[0]  is  x  coordinate,  x[l]  is  y  */ 

RKIND  **z ;  /*  modal  displacements  and  accelerations,  or  displacement 
history  */ 

RKIND  *d; 

RKIND  **res;  /*  residual  */ 

RKIND  *omega;  /*  natural  frequencies,  rad/s  */ 

IKIND  *n;  /*  rows  */ 

IKIND  *m;  /*  columns  */ 

RKIND  *sc;  /*  scale  factor  */ 

RKIND  *h;  /*  time  step  size  */ 

RKIND  *its;  /*  number  of  time  steps  */ 
char  *bl ; 

bl  =  (char  *) chkmalloc (MAXLINE*sizeof (char) ) ; 
sc  =  (RKIND  *) chkmalloc (sizeof (RKIND)) ; 
h  =  (RKIND  *) chkmalloc (sizeof (RKIND) ) ; 
its  =  (RKIND  *) chkmalloc (sizeof (RKIND) ) ; 
n  =  (IKIND  *) chkmalloc (sizeof (IKIND) ) ; 
m  =  (IKIND  *) chkmalloc (sizeof (IKIND) ) ; 
mode_out  =  fopen("mode_hist . dat" , "w") ; 

if  (argc  !  =2)  {printf  ("Usage :  \n  °/0s  input_f  ile_name  \n"  ,  argv  [0] ) ;  exit(l);} 

else{ 

in_file  =  f open(argv  [1] , "r") ; 

if (in_f ile==NULL)  file_error(  argv[l]  ); 

> 

INPUT  *in; 
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in  =  (INPUT  *) chkmalloc (  sizeof (INPUT)  ); 


get_input(  in_file,  in,  &A,  &x,  &z,  &res,  feomega) ; 
d  =  (RKIND  *) chkmalloc (in- >m [0] *sizeof (RKIND) ) ; 
perturb (in,  z,  omega); 
j=0;r=0; 

for(j=0;  j<in->n[0] ; j++){ 
printf ("%g  " , omega [j] ) ; 

> 

printf ("\n") ; 
j=0; 

while(  j  <  in->its [0]  ){/*  time  loop  */ 

/*  write  out  the  shapes  only  every  report  iterations  */ 
if  (  j  %  in->report  [0]  ==  0)-( 

sprintf(bl,  "response-"IF0RM" .out" ,  r) ; 
response  =  fopen(bl,  "w"); 
write_shapes (in,  response,  A,  x,  z,  d) ; 
f close (response) ; 
r++; 

} 

/*  write  the  modal  amplitudes  every  iteration  */ 
write_modes (in,  mode_out,  z) ; 

/*  integrate  */ 
if (  in->method [0]  ==  0  ) 
rk(  in,  z,  omega,  res  ); 
if (  in->method  [0]  ==  1  ) 

semi_implicit (  in,  z,  omega  ); 
if (  in->method [0]  ==  2  ) 

third_order(  in,  z,  omega  ); 
if (  in->method  [0]  ==  3  ) 

fourth_order (  in,  z,  omega  ); 
if (  in->method [0]  ==  4  ) 

fifth_order(  in,  z,  omega  ); 
if (  in->method [0]  ==  5  ) 

sixth_order(  in,  z,  omega  ); 
if (  in->method [0]  ==  6  ) 

eighth_order (  in,  z,  omega  ); 
if (  in->method [0]  ==  7  ) 

coupled_imp_2 (  in,  z,  omega  ); 
if (  in->method [0]  ==  8  ) 

coupled_imp_3 (  in,  z,  omega  ); 
if (  in->method [0]  ==  9  ) 

coupled_imp_4(  in,  z,  omega  ); 
if (  in->method [0]  ==  10  ) 
newmark (  in ,  z ,  omega  ) ; 
if (  in->method [0]  ==  11  ) 
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coupled_var_accel_2(  in,  z,  omega  ); 
if (  in->method  [0]  ==  11  ) 
rk_3(  in,  z,  omega,  res  ); 


j++; 

> 

return(O) ; 


void  write_shapes (INPUT  *in,  FILE  *f,  RKIND  **A,  RKIND  **x,  RKIND  **z,  RKIND  *d) 

{ 

IKIND  i , j ; 

for(j=0;  j  <  in->m[0]  ;  j++){ 
d  [j]  =  0; 

for(i=0; i<in->n[0] ;i++){ 

d [j]  =  d[j]  +  A[i]  [j]*z[0]  [i]  ; 

} 

fprintfCf,  RFORM"  "RFORM"  "RFORM"  \n" ,  x  [0]  [ j ] ,  x [1]  [ j ]  ,  d  [ j ] ) ; 

> 

} 

void  write_modes (INPUT  *in ,  FILE  *f ,  RKIND  **z) 

{ 

IKIND  i; 

for(i=0;  i  <  in->n[0] ;  i++)  fprintf(f,  RFORM"  ",z[0][i]); 
fprintf (f , "\n") ; 


C.5  solve,  c 

#include  " solve. h" 

void  perturb (INPUT  *in ,  RKIND  **z,  RKIND  *omega) 

{ 

IKIND  i; 

for(i=0;  i  <  in->n[0] ;  i++){ 
z[0]  [i]  =  in->perturb  [i] ; 
if (  in->method  [0]  ==  0  ){/*  rk  */ 
z[l]  [i]  =  ZERO; 

}else{ 

if (  in->method  [0]  ==  7  ){/*  coupled  implicit  */ 
z  [1]  [i]  =  in->perturb[i] ; 
z [2]  [i]  =  ZERO; 
z [3]  [i]  =  ZERO; 

}else{ 

if(  in->method [0]  ==  10  ){/*  newmark  */ 
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z[l][i]  =  in->perturb  [i]  ; 
z [2]  [i]  =  ZERO; 
z [3]  [i]  =  ZERO; 
z [4]  [i]  =  ZERO; 
z [5]  [i]  =  ZERO; 

}else-( 

if (in->method[0]  ==  12  ){/*  rk  scheme  with  three  eqns  */ 
z [2]  [i]  =  ZERO; 
z [4] [i]  =  ZERO; 

}else{/*  one  of  the  finite  difference  methods  */ 
z[l] [i]  =  in->perturb [i] ; 
z[2]  [i]  =  in->perturb [i] ; 
z[3][i]  =  in->perturb [i] ; 
z[4][i]  =  in->perturb [i]  ; 
z[5]  [i]  =  in->perturb [i] ; 
z[6]  [i]  =  in->perturb [i] ; 
z[7][i]  =  in->perturb [i] ; 
z[8][i]  =  in->perturb [i] ; 
z[9]  [i]  =  in->perturb [i]  ; 

> 

> 

> 

} 

> 


void  rhs (INPUT  *in,  RKIND  **z,  RKIND  *omega,  RKIND  **res) 
{/*  z[0]  is  displacements,  z[l]  is  velocitiess  */ 

IKIND  j ; 

f or (j=0;  j  <  in->n[0]  ;  j++){ 
res  [0]  [j]  =  z[l]  [j]  ; 
res  [1]  [j]  =  -omega [j]  *z[0]  [j]  ; 

> 

> 

void  rk(INPUT  *in,  RKIND  **z,  RKIND  *omega,  RKIND  **res) 

{ 

IKIND  i , j , s ; 

RKIND  * alpha; 

IKIND  stages  =  4; 

RKIND  **zl ; 

zl  =  (RKIND  **) chkmalloc (2*sizeof (RKIND) ) ; 
zl [0]  =  (RKIND  *) chkmalloc (in->n [0] *sizeof (RKIND) ) ; 
zl [1]  =  (RKIND  *) chkmalloc (in->n [0] *sizeof (RKIND) ) ; 
alpha  =  (RKIND  *) chkmalloc (stages*sizeof (RKIND) ) ; 
alpha [0]  =  ONE/FOUR; 
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alpha [1]  =  ONE/THREE; 
alpha [2]  =  ONE/TWO; 
alpha [3]  =  ONE; 

/*  store  the  old  data  */ 
f or (j=0;  j  <  2;  j++H 

for(i=0;  i  <  in->n[0] ;  i++){ 
zl[j]  [i]  =  z[j]  [i]  ; 

} 

> 

for(s=0;  s<stages;  s++){ 

rhs(  in,  z,  omega,  res); 

f or (j=0;  j  <  2;  j++){ 

for(i=0;  i  <  in->n[0] ;  i++){ 

z[j][i]  =  zl  [j]  [i]  -  alpha  [s]  *res  [j]  [i]  *in->h[0]  ; 

> 

} 

> 

> 

void  newmarkC INPUT  *in,  RKIND  **z,  RKIND  Iomega) 

{/*  z:  0,  1  displacement 
2,  3  velocity 
4,  5  acceleration  */ 

RKIND  dt; 

IKIND  i; 

dt  =  in->h[0]  ; 

for(  i=0;  i<in->n[0] ;  i++  ){ 

/*  slide  the  solution  back  */ 
z  [5]  [i]  =  z  [4]  [i]  ; 
z  [3]  [i]  =  z  [2]  [i]  ; 
z[l]  [i]  =  z  [0]  [i]  ; 

/*  update  acceleration  */ 

z[4][i]  =  -omega  [i]  *  (  z  [0]  [i]  +  dt*z[2][i]  +  dt*dt*z  [5]  [i]  /3  )/ 

(  1  +  dt*dt*omega[i] /6  ); 

/*  update  velocity  */ 

z [2]  [i]  =  z [3]  [i]  +  dt*  (  z [5]  [i]  +  z [4]  [i]  )/TW0; 

/*  update  displacement  */ 

z[0]  [i]  =  z[l]  [i]  +  dt*z[3]  [i]  +  dt*dt*(  TW0*z[5]  [i]  +  z[4]  [i]  )/SIX 

> 

} 
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void  coupled_var_accel_2 (INPUT  *in ,  RKIND  **z,  RKIND  Iomega) 

{/*z:  0,1  displacement 
2,3  velocity 
4,5  acceleration  */ 

RKIND  a,  bl,  b2,  b3,  dt; 

IKIND  i; 

dt  =  in->h[0]  ; 

for(i=0;  i<in->n[0] ;  i++){ 

a  =  FOUR*dt*dt*omega[i]  +  9.0; 

/*  right  hand  side  vector  */ 
bl  =  F0UR*z  [0]  [i]  -  z[l]  [i]  ; 

b2  =  F0UR*z  [2]  [i]  -  z  [3]  [i]  ; 

b3  =  F0UR*z  [4]  [i]  -  z  [5]  [i]  ; 

/*  move  the  solution  back  */ 

z[l]  [i]  =  z  [0]  [i]  ; 
z  [3]  [i]  =  z  [2]  [i]  ; 
z  [5]  [i]  =  z  [4]  [i]  ; 

/*  calculate  the  new  solution  */ 

z [0]  [i]  =  bl /THREE  +  TW0*dt*b2/a  +  F0UR*dt*dt*b3/ (3*a) ; 

z  [2]  [i]  =  THREE*b2/a  +  TW0*dt*b3/a; 

z[4][i]  =  -TW0*dt*omega[i] *b2/a  +  THREE*b3/a; 

> 

} 

void  semi_implicit (INPUT  *in,  RKIND  **z,  RKIND  *omega) 

{ 

IKIND  i; 

/*  move  everything  back  one  */ 
for(i=0;  i<in->n[0] ;  i++){ 
z  [3]  [i]  =  z  [2]  [i]  ; 
z  [2]  [i]  =  z[l]  [i]  ; 
z[l]  [i]  =  z  [0]  [i]  ; 

> 

for(  i=0;  i<  in->n[0] ;  i++  ){/*  loop  over  the  modes  */ 

z[0][i]  =  (  FIVE*z[l][i]  -  F0UR*z[2][i]  +  z [3]  [i]  )/(  TWO  +  pow(in->h[0]  , TWO) *omega [i] 

> 

} 

void  third_order (INPUT  *in,  RKIND  **z,  RKIND  *omega) 

{/*XXX  unstable  XXX*/ 

IKIND  i; 

RKIND  a; 

/*  move  everything  back  one  */ 
for(i=0;  i<in->n[0] ;  i++){ 
z  [4]  [i]  =  z  [3]  [i]  ; 
z  [3]  [i]  =  z  [2]  [i]  ; 
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z[2]  [i]  =  z[l]  [i]  ; 
z[l]  [i]  =  z  [0]  [i]  ; 

> 

for(  i=0;  i<  in->n[0] ;  i++  ){/*  loop  over  the  modes  */ 
a  =  (  35.0  +  12 . 0*pow(in->h [0] , TWO) Iomega [i]  ); 

z  [0]  [i]  =  (  104 . 0*z  [1]  [i]  -  114 . 0*z  [2]  [i]  +  56 . 0*z  [3]  [i]  -11 . 0*z  [4]  [i]  ) 

/  a; 

> 

> 

void  fourth_order( INPUT  *in,  RKIND  **z,  RKIND  Iomega) 

{/*XXX  unstable  XXX*/ 

IKIND  i; 

RKIND  a; 

/*  move  everything  back  one  */ 
for(i=0;  i<in->n[0] ;  i++){ 
z  [5]  [i]  =  z  [4]  [i]  ; 
z  [4]  [i]  =  z  [3]  [i]  ; 
z  [3]  [i]  =  z  [2]  [i]  ; 
z  [2]  [i]  =  z[l]  [i]  ; 
z[l]  [i]  =  z  [0]  [i]  ; 

> 

for(  i=0;  i<  in->n[0] ;  i++  ){/*  loop  over  the  modes  */ 
a  =  (  45.0  +  12 . 0*pow(in->h [0] ,TW0) *omega[i]  ); 

z  [0]  [i]  =  (  154 . 0*z  [1]  [i]  -  214 . 0*z  [2]  [i]  +  156 . 0*z  [3]  [i]  -  61 . 0*z  [4]  [i]  +  10 . 0*z  [5]  [i] 

> 

} 


void  fifth_order (INPUT  *in,  RKIND  **z,  RKIND  *omega) 
{/*XXX  this  one  is  unstable  XXX*/ 


IKIND  i; 
RKIND  a; 


/*  move  everything  back  one  */ 
for(i=0;  i<in->n[0] ;  i++){ 


z [6]  [i] 

=  z  [5] 

[i] ; 

z [5]  [i] 

=  z  [4] 

[i] ; 

z [4]  [i] 

=  z  [3] 

[i] ; 

z [3]  [i] 

=  z  [2] 

[i] ; 

z [2]  [i] 

=  z  [1] 

[i] ; 

z[l]  [i] 

=  z  [0] 

[i] ; 

> 

for(  i=0; 

i<  in- 

>n[0] ;  i++  ){/*  loop  over  the  modes  */ 

a  =  (  203/45  +  pow(in->h[0] ,TW0)*omega[i]  ); 

z  [0]  [i]  =  (  87*z[l][i]/5  -  117*z  [2]  [1] /4  +  254*z  [3]  [i] /9  -  33*z[4][i]/2 
+  27*z[5][i]/5  -  137*z  [6]  [i]/180  )/a; 
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> 

> 

void  sixth_order (INPUT  *in,  RKIND  **z,  RKIND  *omega) 

{ 

IKIND  i; 

RKIND  a; 

/*  move  everything  back  one  */ 
for(i=0;  i<in->n[0] ;  i++){ 
z[7]  [i]  =  z  [6]  [i]  ; 
z[6]  [i]  =  z  [5]  [i]  ; 
z[5]  [i]  =  z  [4]  [i]  ; 
z  [4]  [i]  =  z  [3]  [i]  ; 
z  [3]  [i]  =  z  [2]  [i]  ; 
z  [2]  [i]  =  z[l]  [i]  ; 
z[l]  [i]  =  z  [0]  [i]  ; 

> 

for(  i=0;  i<  in->n[0] ;  i++  ){/*  loop  over  the  modes  */ 
a  =  (  938  +  180*pow(in->h [0] , TWO) Iomega [i]  ); 

z [0]  [i]  =  (  4014*z  [1]  [i]  -  7911*z  [2]  [i]  +  9490*z  [3]  [i]  -  7380*z[4][i] 
+  3618*z  [5]  [i]  -  1019*z  [6]  [i]  +  126*z[7][i]  )/a; 

> 

> 

void  eighth_order (INPUT  *in,  RKIND  **z,  RKIND  Iomega) 

{/*XXX  unstable  XXX*/ 

IKIND  i; 

RKIND  a; 

/*  move  everything  back  one  */ 
for(i=0;  i<in->n[0] ;  i++){ 
z  [9]  [i]  =  z  [8]  [i]  ; 
z  [8]  [i]  =  z  [7]  [i]  ; 
z  [7]  [i]  =  z  [6]  [i]  ; 
z  [6]  [i]  =  z  [5]  [i]  ; 
z  [5]  [i]  =  z  [4]  [i]  ; 
z  [4]  [i]  =  z  [3]  [i]  ; 
z  [3]  [i]  =  z  [2]  [i]  ; 
z  [2]  [i]  =  z[l]  [i]  ; 
z[l]  [i]  =  z  [0]  [i]  ; 

> 

for(  i=0;  i<  in->n[0] ;  i++  ){/*  loop  over  the  modes  */ 
a  =  (  32575  +  5040*pow(in->h [0] ,TW0)*omega[i]  ); 
z [0]  [i]  =  (  165924*z  [1]  [i]  -  422568*z  [2]  [i]  +  704368*z  [3]  [i] 

-  818874*z  [4]  [i]  +  667800*z  [5]  [i]  -  375704*z  [6]  [i] 

+  139248*z  [7]  [i]  -  30223*z  [8]  [i]  +  3044*z[9][i]  )/a; 

> 
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> 


void  coupled_imp_2 (INPUT  *in ,  RKIND  **z,  RKIND  Iomega) 
{/*z  0,  1  :  displacement 
2,  3  :  velocity  */ 

IKIND  i; 

RKIND  a,  bl,  b2; 

for(i=0;  i  <in->n[0] ;  i++){/*  loop  over  the  modes  */ 
a  =  F0UR*in->h [0] *in->h[0] Iomega [i]  +  9; 
bl  =  F0UR*z  [0]  [i]  -  z[l]  [i]  ; 
b2  =  F0UR*z  [2]  [i]  -  z  [3]  [i]  ; 

/*  move  the  velocity  back  one  */ 
z  [3]  [i]  =  z  [2]  [i]  ; 

/*  modve  the  displacements  back  one  */ 
z[l]  [i]  =  z  [0]  [i]  ; 

/*  calculate  the  new  ones  */ 

z [0]  [i]  =  (  THREE*bl  +  2*in->h [0] *b2  )/a; 

z [2]  [i]  =  (  -TW0*in->h [0] *omega [i] *bl  +  THREE*b2  )/a; 

> 

> 

void  coupled_imp_3 (INPUT  *in,  RKIND  **z,  RKIND  Iomega) 
{/*XXX  unstable  XXX*/ 

/*z  0,  1,  2  :  displacement 
3,  4,  5  :  velocity  */ 

IKIND  i; 

RKIND  a,  bl,  b2; 

for(i=0;  i  <in->n[0] ;  i++){/*  loop  over  the  modes  */ 
a  =  36 . 0*in->h  [0] *in->h[0] *omega [i]  +  121.0; 
bl  =  18 . 0*z  [0]  [i]  -  9 . 0*z  [1]  [i]  +  2 . 0*z  [2]  [i]  ; 
b2  =  18 . 0*z  [3]  [i]  -  9 . 0*z  [4]  [i]  +  2. 0*z  [5]  [i]  ; 

/*  move  the  displacements  back  one  */ 
z  [2]  [i]  =  z[l]  [i]  ; 
z[l]  [i]  =  z  [0]  [i]  ; 

/*  move  the  velocities  back  one  */ 
z  [5]  [i]  =  z  [4]  [i]  ; 
z  [4]  [i]  =  z  [3]  [i]  ; 

/*  calculate  the  new  ones  */ 

z [0]  [i]  =  (11 . 0*bl  +  6 . 0*in->h[0] *b2) /a; 

z[3][i]  =  (-6 . 0*in->h  [0] *omega[i] *bl  +  11.0*b2)/a; 

> 

> 

void  coupled_imp_4 (INPUT  *in,  RKIND  **z,  RKIND  *omega) 
{/*XXX  unstable  XXX*/ 

/*z  0,  1,  2,  3  :  displacement 
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4,  5,  6,  7  :  velocity  */ 

IKIND  i; 

RKIND  a,  bl,  b2; 

for(i=0;  i  <in->n[0] ;  i++){/*  loop  over  the  modes  */ 
a  =  144*in->h[0] *in->h[0] Iomega [i]  +  625; 
bl  =  48*z  [0]  [i]  -  36*z[l][i]  +  16*z[2][i]  -3*z[3][i]; 
b2  =  48*z  [4]  [i]  -  36*z  [5]  [i]  +  16*z[6]  [i]  -3*z[7]  [i]  ; 
/*  move  the  displacements  back  one  */ 
z  [3]  [i]  =  z  [2]  [i]  ; 
z  [2]  [i]  =  z[l]  [i]  ; 
z[l]  [i]  =  z  [0]  [i]  ; 

/*  move  the  velocities  back  one  */ 
z  [7]  [i]  =  z  [6]  [i]  ; 
z  [6]  [i]  =  z  [5]  [i]  ; 
z  [5]  [i]  =  z  [4]  [i]  ; 

/*  calculate  the  new  ones  */ 

z [0]  [i]  =  (25*bl  +  12*in->h[0] *b2)/a; 

z[4][i]  =  (-12*in->h  [0] *omega [i] *bl  +  25*b2)/a; 

> 

> 

void  rk_3(INPUT  *in,  RKIND  **z,  RKIND  *omega,  RKIND  **res) 
{/*z:  0,1  displacement 
2,3  velocity 
4,5  acceleration  */ 

IKIND  i , j ; 

RKIND  dt; 

RKIND  * alpha; 

IKIND  stages  =  4; 

alpha  =  (RKIND  *) chkmalloc (stages*sizeof (RKIND) ) ; 

dt  =  in->h[0]  ; 
alpha [0]  =  ONE/FOUR; 
alpha [1]  =  ONE/THREE; 
alpha [2]  =  ONE/TWO; 
alpha [3]  =  ONE; 

/*  store  the  old  solution  */ 
for(i=0; i<in->n[0] ;i++){ 
z[l]  [i]  =  z  [0]  [i]  ; 
z  [3]  [i]  =  z  [2]  [i]  ; 
z  [5]  [i]  =  z  [4]  [i]  ; 

> 

for(j=0;  j  <  stages;  j++){ 
for(i=0;  i  <  in->n[0] ;  i++){ 
res  [0]  [i]  =  z  [2]  [i]  ; 
res  [1]  [i]  =  z  [4]  [i]  ; 
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res [2]  [i]  =  -omega  [i] *z  [2]  [i]  ; 
z[0][i]  =  z[l][i]  -  alpha[j] *dt*res  [0]  [i] 
z[2][i]  =  z[3][i]  -  alpha[j] *dt*res [1] [i] 
z[4][i]  =  z[5][i]  -  alpha[j] *dt*res [2] [i] 

} 

> 

> 
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Appendix  D.  Fast  Fourier  Transform 

This  routine  uses  the  fftw  library  released  by  MIT  under  the  GNU  General  Public  License, 
and  available  at  http://www.fftw.org. 

D.  1  cl-jft.  c 

/*  This  is  a  simple  program  that  reads  a  column  of  data  from  a  text 

file  and  performs  a  discrete  fourier  transform  on  it,  normalizes  the 
output  by  the  DC  offset,  (if  it  isn't  zero),  outputs  the  result  to 
another  text  file  */ 

#include  <stdlib.h> 

#include  <stdio.h> 

#include  <math.h> 

#include  <fftw3.h> 

#include  "getwords.h" 

#def ine  MAXLINE  200 
#def ine  MAXWORDS  2 

int  main(  int  argc,  char  *argv[]) 

{ 

double  *f,  dt,  norm; 
double  *in; 
fftw_complex  *out; 
fftw_plan  p; 
int  n,  ac,  i; 

FILE  *in_f ; 

FILE  *out_f; 
char  line [MAXLINE] ; 
char  *av [MAXWORDS] ; 
if (  argc  !=  3  ){ 

printf ( "Error . \nUsage : \n" 

"  Is  SAMPLE_RATE  dataf ile\n" , argv [0] ) ; 
exit (1) ; 

> 

dt  =  atof (  argv[l]  ); 

in_f  =  fopen(  argv [2],  "r"  ); 

sprintf (  line,  "70s.fft",  argv [2]  ); 

out_f  =  fopen(line,  "w"); 

if (  in_f  ==  NULL  ){ 

printf  ("Error\nCould  not  open\n  °/0s\nfor  read.  \n" , argv [1] ) ; 
exit (1) ; 

> 

n  =  0; 

while (  f gets (line,  MAXLINE,  in_f  )  !=  NULL  ){ 
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n++; 

> 

if (  n  ==  0  ){ 

printf ("Error . \nNo  lines  of  data.\n"); 
exit (1) ; 

> 

fclose(in_f) ; 

in  =  (double  *)malloc(  sizeof (double) *n  ); 
f  =  (double  *)malloc(  sizeof (double) *n/2  ); 

out  =  (f ftw_complex  *)fftw_malloc(  sizeof (fftw_complex) * (l+n/2)  ); 
in_f  =  fopen(  argv[2],  "r") ; 
n  =  0; 

while (  f gets (line,  MAXLINE,  in_f  )  !=  NULL  ){ 
ac  =  getwords (line ,  av,  MAXW0RDS-1  ); 
in[n]  =  atof  (  av[0]  ); 
n++; 

> 

p  =  f f tw_plan_dft_r2c_ld(  n,  in,  out,  FFTW_ESTIMATE  ); 
f f tw_execute (p) ; 
for(i=0; i<n/2; i++){ 

/*  frequency  real  imaginary  magnitude  */ 
if(i==0){ 

if  (  out  [i]  [0]  !  =  0  ) 
norm  =  out  [i]  [0] ; 

else 
norm  =  1.0; 

printf  ("Average  :  7og\n" ,  norm/n) ; 

> 

in[i]  =  sqrt(  out [i] [0] *out  [i]  [0]  +  out [i] [1] *out  [i]  [1]  ); 
f [i]  =  ( (double) i) / (  dt*(double)n  ); 

f printf  (out_f ,  "°/0g  °/,g  %g  %g\n"  ,  f  [i]  ,  out  [i]  [0] /norm,  out  [i]  [1] /norm, 
in [i] /norm  ); 

> 

f f tw_destroy_plan(p) ; 
f f tw_f ree (in) ;  f f tw_free (out) ; 

return(0) ; 
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And  now  for  something  completely  different  ... 


Sonnet  of  a  MAV 


We  know  flapping  wings,  and  flopping  wings. 

That  ought  drive  us  to  distraction. 

Isn’t  it  always  the  membrane  things. 

That  create  the  main  attraction? 

Their  carbon  fiber  ribs  and  spars 
With  nylon  in  betwixt 
Refuse  to  bend  or  stretch  too  far 
With  fluids  in  the  mix. 

And  towards  the  middle  of  the  night. 

When  most  folks  are  gone  to  sleep, 

We  let  our  bits  in  machines  take  flight, 

So  we  know  those  to  throw  away  or  keep. 

But  in  the  last  analysis  our  derivations  were  in  vain. 
The  world,  it  is  non-linear!  Now  isn’t  that  a  pain? 
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